1  .  RD-A149  956  OPTIMIZATION  OF  GUIDANCE  AND  CONTROL 

1  MINIMIZATION  AND  NAVSTAR/GPSCU)  NAVAL 
1  SCHOOL  MONTEREV  CA  V  C  GARCIA  SEP  84 
1  UNCLASSIFIED 

USING  FUNCTION 
POSTGRADUATE 

F/G  20/2 

1/1 

NL 

^  k 

_ 1 

NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


33B^3CT2HSB!EBZIiniIIi 


REPORT  DOCUMENTATION  PAGE 


4.  TITLE  (m<t  SubUtI*) 


Optimization  of  Guidance  and  Control  using 
Function  Minimization  and  NAVSTAR/GPS 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


RECIPIENT'S  CAT ALOO  NUMBER 


S.  TYRE  OP  REPORT  4  PERIOD  COVERED 


Master's  Thesis 
September  1984 


«.  PERPORMINO  ORO.  REPORT  NUMBER 


7.  AUTHORS 

Vicente  Chavez  Garcia,  Jr. 


PERPORMINO  ORGANISATION  NAME  AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California  93943 


II.  CONTROLLING  OPPICE  NAME  ANO  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California  93943 


.  MONITORING  AGENCY  NAME  A  ADDRESSfff  dlllotoni  I 


OMTAACT  OR  GRANT  NUMGERfaJ 


I 


12.  REPORT  OATE 

September  1984 _ 


IS.  NUMBER  OP  PAGES 

94 


i  Controlling  Otllea)  IS.  SECURITY  CLASS,  (el  thlo  report; 

UNCLASSIFIED  _ 

Mo.  OECLASSIPICATION/OOWNGRADINO 
SCHEDULE 


IS.  DISTRIBUTION  STATEMENT  (ol  thlo  Report; 


Approved  for  public  release,  distribution  unlimited 


19.  KEY  WORDS  (Continue  on  rovotoo  ml  do  II  noeoooory  end  Identity  by  block  number) 


Optimization,  Performance  Criterion,  Narroto,  Taylor's  Series, 

Function  Minimization  Subroutine,  NAVSTAR/GPS,  Space  Transportation  System 
(STS), . Sea-Land  Mclean  (SL-7)^  / 

L  J 


0.  ABSTRACT  (Conllnuo  on  reverie  tldo  If  nocoooorr  and  Identity  by  black  nuetbor) 

A  carefully  designed  controller,  tuned  to  minimize  a  performance  criterion 
based  on  representation  of  the  added  drag  due  to  steering,  can  minimize 
propulsion  losses.  A  computer  simulation  modeling  the  Sea-Land  Mclean  (SL-7) 
oontauLnership  was  coupled  to  a  function  minimization  subroutine  and  a  sea- 
state  generator  subroutine  to  accomplish  the  tuning.  Storing  these  optimal 
controller  parameters  in  a  look  up  table  as  functions  of  ship  state,  sea 
state,  and  encounter  angle,  this  technique  can  be  used  as  an  adaptive 


I  JAN  72 


COITION  OP  I  NOV  ••  12  OBSOLETE 
S/N  0102-  LP-  014-  6401 


SECURITY  CLASSIFICATION  OP  THIS  PABE  | 


JJisr 


20.  Continued 

controller.  Satellite  platforms  can  give  continuous  environmental  operating 
conditions  which  may  be  used  to  select  proper  controller  parameters  to 
provide  continuous  operation  on  a  mininun  of  the  cost  function.  The  SL-7 
contai  nership  on  a  minimum  of  the  cost  function.  The  SL-7  containership 
computer  model  was  tested  in  calm  waters  and  in  a  seaway. 


nt  y  i 

\,At  C  L  u/O  £  ' 


't  f  £3 


>  f  y 


SfCURITV  CLASSIFICATION  OF  THIS 


S  N  0102-  LF*  014*  6601 


Approved  for  public  release;  distribution  unlimited 


Optiaiz?tion  of  $uid?nce  and  Control  using 
Function  Bimmzation  ana  NAVSTAB/GPS 

by 


Vicente  Chavez  Garcia,  Jr. 
Lieutenant,  United  States  Navy 
B.S.E.E.,  New  Mexico  State  University,  1978 
B.S.E. ,  University  of  Central  Florida,  1982 


Subnitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 

fron  the 

NAVAL  POSTGRADUATE  SCHOOL 
Septeaber  1984 


I ^2rltd '  pv5/?^W  -ffru  //•  R 

Harriett ~B7-Hi?5gT~cEairisB7LDgpjrtfl6Bf~ 

Electrical  and  Coaputer  Engineering 

- ;"^/^^9^n-HrDyef7“". - - - 

Eean  of  Science  ana  Engineering 


A  carefully  designed  controller,  tuned  to  minimize  a 
performance  criterion  based  on  representation  of  the  added 
drag  due  to  steering,  can  nininize  propulsion  losses.  A 
coaputer  simulation  modeling  the  Sea-Land  Hclean  (SL-7) 
containership  vas  coupled  to  a  function  minimization  subrou¬ 
tine  and  a  sea-state  generator  subroutine  to  accomplish  the 
tuning.  Storing  these  optimal  controller  parameters  in  a 
look  up  table  as  functions  of  ship  state,  sea  state,  and 
encounter  angle,  this  technique  can  be  used  as  an  adaptive 
controller.  Satellite  platforms  can  give  continuous  environ¬ 
mental  operating  corditions  which  may  be  used  to  select 
proper  controller  parameters  to  provide  continuous  operation 
on  a  minimum  of  the  cost  function.  The  SL-7  containership 
computer  model  vas  tested  in  calm  waters  and  in  a  seaway. 


TABLE  OF  CONTENTS 


I.  I NTBODUCT ION . 

II.  COMPUTES  SODELS  . 

III.  COST  FUNCTION  ..  . 

IV.  CCNT80LLEH  DESIGN  USING  ICSOS  . 

V.  CCNTBOLLEB  DESIGN  USING  FOBTBAN  PBOGB&N  . 

VI.  CCNTBOLLEB  DESIGN  IN  SEA  STATE  . 

VII.  AN  ADAPTIVE  CCNTBOLLEB . ‘  .  . 

VIII.  CONCLUSIONS  1ND  BECOHHENDATIONS  FOB  FUTUBE 

STUDY  . 

APPENDII  A;  PBOGBAH  TO  CALCULATE  OPTIMAL  GAINS  . 

APPENDIX  E:  EXAMPLE  PBOBLEM  USING  ICSOS  .  .  .  . 

LIST  OF  BEFEBENCES  . 

INITIAL  DISTRIBUTION  IIST  . 


i  Accession  For 

NTTS  -1 
rii"  t  ” 


LIST  OF  TABLES 


THIBD-OHDEB  NCHOTO 

HODEL  FOB  THE  SL-7 

•  •  • 

• 

• 

9 

15 

HYDBODYNAHIC  COEFFICIENTS  FOB  THE 

SL-7 

•  •  m 

m 

• 

16 

BEIGHTING  FACTCB  . 

• 

•  •  •  • 

•  •  « 

•  • 

•  mm 

• 

m 

• 

18 

BEID'S  BESOLTS  .  . 

m 

•  •  •  • 

•  •  « 

•  • 

•  mm 

9 

m 

• 

20 

ICSOS  BESOLTS  .  . 

• 

•  •  •  • 

•  •  « 

•  • 

•  •  • 

m 

m 

• 

20 

SIflOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

m 

• 

• 

22 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

22 

SIflOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

m 

• 

22 

SIMULATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

23 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

23 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

m 

• 

24 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

9 

m 

• 

25 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

m 

• 

27 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

27 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

m 

9 

• 

28 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

9 

29 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

m 

9 

30 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

42 

SIHOLATION 

RESULTS 

- 

STEADY 

STATE 

600 

SECS 

9 

• 

• 

42 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

m 

• 

• 

43 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

43 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

9 

• 

• 

44 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

• 

• 

44 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

• 

m 

• 

45 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600 

SECS 

m 

m 

• 

45 

SIHOLATION 

BESOLTS 

- 

STEADY 

STATE 

600  SECS  . 

m 

• 

9 

46 

CCHPABISON 

OF  THE  HINIHOH  COST  . 

•  • 

•  •  • 

• 

• 

• 

46 

EFFECTS  DDE  TO  TRANSIENT  AND  GHADDAL  BUILD  OP 

Of  SEA  STATE . 47 


48 

49 
49 


SIMULATION  RESULTS  -  STEAD!  STATE  600  SECS  ...  48 

SIMULATION  BESGLTS  -  STEAD!  STATE  600  SECS  ...  49 

SIMULATION  BESULTS  -  STEAD!  STATE  600  SECS  ...  49 

COHPABISON  OF  IBBEGULAB  TO  BEGULAB  SEAS 

CCNTBOLLEB  GAINS  .  50 

CCMPABISON  OF  IBBEGULAB  TO  BEGULAB  SEAS 

CCNTBOLLEB  GAINS  . .50 

ICSOS  OUTPUT . 90 


r-1 


■  1 

*  .  J 

\*/*j 

•  •  •  *  *  »  ■  *  •  -  » 


LIST  OF  FIGURES 


BLOCK  DIAGHJfl  . 

DETERMINATION  OF  NOHOTO  MODELS  . 

OPTIMIZATION  OF  CONTSOLLEB  . 

VARIOUS  STB  OCT (JRES  FOB  CONTROLLERS  .  .  . 
OPTIMIZATION  OF  STATE  FEEDBACK  CONTROLLER 
OPTIMIZATION  OF  CONTBCLLEB  USING  FORTRAN 

PROGRAM  . 

THO  STATE  SYSTEM . 

THREE  STATE  SYSTEM  . 

YAH  vs.  TIME  (controller  A,  B  j  C  and 
state-feedback)  ............ 

RODDER  vs.  TIME  (controller  A,  B,  C  and 

state- feedback)  . . . 

YAH  AND  RUDDER  vs.  TIME  (controller  A) 

YAH  AND  RODDER  vs.  TIME  (controller  B) 

YAH  AND  RODDER  vs.  TIME  (controller  C) 

YAH  AND  RUDDER  vs.  TIME  (state-feedback 

controller)  .............. 

OPTIMIZATION  OF  CONTROLLER  IN  SEA  STATE 
ADDED  NASS  VS.  ENCOUNTER  FREQUENCY  .  .  . 
ADDED  INERTIA  vs.  ENCOUNTER  FREQUENCY  . 

ENERGY  DENSITY  SPECTRON  . 

XAH  VS.  TIME  30  DEGREES  . 

RUDDER  VS.  TIME  30  DEGREES  . 

YAH  VS.  TIME  60  DEGREES  . 

RUDDER  VS.  TIME  60  DEGREES  . 

IAH  vs.  TINE  90  DEGREES  . 

RUDDER  VS.  TIME  90  DEGREES  . 


TAB  vs.  TI HI  120  DEGREES  . 
fiOOOBB  VS.  TIME  120  DEGREES 
TAB  vs.  TINE  150.  DEGBEES  . 
BODDEB  VS.  TISE  150  DEGBEES 
ADAPTIVE  CON1BOLLEB  .  .  .  . 
DIGITAL  BLOCK  DIAGBAH  .  .  . 
HATCHED  PILTEB  DESIGH  .  .  . 
DETAIL  BLOCK  DIAGBAM  .  .  .  . 
TAB  AND  BUDDEB  VS.  TINE  .  . 


ACKHOB LEDGE HEHT 


r*1 


This  thesis  vas  aade  possible  by  the  continuous  support 
of  ay  wife,  Eailia.  I  thank  iJCDR  Jia  Cass  for  installing  and 
debugging  the  sea  state  generator  prograa  obtained  fron 
DTHSHDC.  I  also  give  thanks  to  LT  Eaaanuel  Horianopoulos  and 
LT  Pericles  Kyritsis  Spyroailios  for  their  assistance  in 
generating  plots  and  data.  Finally,  I  owe  a  great  debt  of 
gratitude  to  Dr.  George  J.  Thaler,  for  his  support,  enthu¬ 
siast,  and  outstanding  exanple  of  excellence. 


I.  IHTBODOCTION 


An  overall  rise  in  fuel  prices  has  led  to  an  increasing 
interest  in  the  design  of  autopilots  for  ships.  The  purpose 
of  the  automatic  steering  control  is  to  minimize  the  propul¬ 
sion  losses,  vhich  are  caused  by  added  drag  due  to  steering 
of  the  ship.  Minimizing  a  performance  criterion  based  on 
added  drag  due  to  steering  can  reduce  fuel  consumption. 
Claims  by  many  researchers  indicate  that  a  carefully 
designed  controller  could  save  from  one  to  two  percent  of 
fuel.  For  large  containerships  this  could  amount  to  more 
than  $100,000.00  per  year  savings. 

To  study  the  optimization  problem,  models  of  both  the 
ship  and  its  operating  environment  are  required.  Shat  type 
of  computer  model  should  be  used  to  represent  the  ship? 
Chapter  two  addresses  the  development  of  several  models. 
Since  the  best  model  mas  desired  it  was  decided  to  use  the 
equations  of  motion  to  simulate  the  ship  in  our  Fortran 
program.  The  basic  Nomoto  models  give  an  adequate  descrip¬ 
tion  of  ship  steering  dynamics  for  design.  The  Ncmoto 
second-  and  third-order  models  were  developed  from  the  equa¬ 
tions  of  motion  as  defined  by  a  series  expansion  including 
all  terms  (both  linear  and  nonlinear)  for  which  hydrodynamic 
coefficients  were  available.  An  interactive  program  that 
utilize'!  the  Nomoto  models  to  model  the  ship  was  also  used. 
Two  independent  programs  were  developed  to  aid  in  the  design 
of  the  controller. 

Bhat  is  an  adequate  cost  function  which  represents  the 
added  drag  due  to  steering?  Chapter  three  addresses  the 
classical  cost  function  used  by  many  researchers. 

Since  a  variety  of  control  algorithms  are  possible  one 
must  ash  if  one  algorithm  provides  a  lower  minimal  cost  than 
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another.  Chapters  four,  five  and  six  address  the  selection 
of  the  controller  which  provides  the  minimum  value  of  added 
drag  due  to  steering. 

Ship  dynamics  change  with  operating  conditions  such  as 
ship  speed  ,  sea  state  ,  and  encounter  angle.  Therefore  an 
adaptive  controller  oust  be  used  to  provide  minimum  added 
drag  due  to  steering.  Chapter  seven  development  of  an 
approach  to  an  adaptive  controller  utilizing  satellite 
information. 

Conclusions  were  drawn  from  these  experiments  and  are 
presented  in  Chapter  eight.  This  thesis  investigated  only 
course  keeping  with  emphasis  on  minimizing  rudder  and  yawing 
activity  to  reduce  fuel  consumption.  Presented  in  this 
Chapter  are  recommendations  for  future  study  where  the 
objective  is  track  following  which  would  be  important  for 
ships  reguired  to  follow  stringent  routes.  It  is  also  impor¬ 
tant  for  other  systems  such  as  satellites,  missiles, 
aircraft,  where  the  controller  minimizes  yaw  error  to  keep 
the  system  on  track. 
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II,  CO HP OT IB  HODELS 

The  aodel  vhich  test  represents  ship-steering  djnaeics 
is  a  Taylor’s  series  expansion  of  the  force  and  aoaent  rela¬ 
tionships  around  a  selected  steady-state  operating  point. 
The  resulting  eguations  are  coaaonly  known  as  the  equations 
of  notion  £Bef.  1  ].  A  coaputer  program  was  developed  using 
known  available  data  on  the  hydrodynaaic  coefficients  for 
the  SL-7  containership  to  provide  a  computer  simulation  of 
the  ship.  The  coaputer  program  is  shown  in  Appendix  A. 
figure  2,1  shows  the  block  diagram.  Small  yaw  command 
angles  are  used,  for  exaaple  7AHC-  1.0  /  57.296  represents  a 
yaw  ccamand  change  of  one  degree. 
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figure  2.1  BLOCK  DI1GBAH 
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To  obtain  the  Moioto  second-  and  third-order  transfer 
functions  fron  the  equations  of  aotionr  the  function  aini- 
nization  subroutine  was  used  to  obtain  the  coefficients, 
figure  2.2  shovs  the  scheae  used  to  obtain  the  Hcaoto 
transfer  functions.  The  computer  prograa  is  shown  in 
Appendix  A. 


figure  2.2  DSTZ1HXIATI0B  Of  IOBOXO  HODELS 

The  lonoto  aodels  were  checked  against  analytic  results 
froa  linearized  equations. 

Proceeding  to  the  second-order  Hoaoto  equation: 


1« 


*(S)/a(S)  =  k/s*  (i  ♦i*s) 


(2.1) 


Deriving  the  second-ccder  Hoaoto  transfer  function  froa  the 
yan  equation  only,  the  resalt  is 

*(S)/«{S)  =  0. 040893/S* ( 1+8. 539932*S) 
and  using  function  ainiaization  as  in  Figure  2.2 

*(S)/«(S)  =  0.040922 1/S*  (1*8. 5520782*S) 
and  the  agreeaent  is  obvious.  Osing  function  ainiaization 
with  both  yaw  and  sway  equations  with  linear  teras  only,  the 
results  are: 

*(S)/*(S)  =  0.1072741/S*  (1*31. 9199524*S) 

If  the  nonlinear  teras  are  included  hut  the  perturbation  is 
saall 

*(S)/a(S)  a  0.  1072082/S*  (1*31. 8907013*S) 
and  it  is  clear  that  the  nonlinear  terns  contribute  little. 

Proceeding  to  the  third- order  Nonoto  equation: 

*(S)/*(S)  =  K*{UTZ*S)/S*  (1*TP1*S)  *  (1*TP2*S)  (2.2) 

The  paraaeters  were  calculated  and  checked  by  using  function 
ainiaization  as  in  Figure  2.2.  The  results  are  given  in 
Table  1.  It  is  clear  that  the  answers  obtained  by  function 
ainiaization  agree  closely  with  the  analytic  solutions. 


TABLE  1 

TEIBD-OBDEB  IOHOTO  BODEL  FOB  THE  SL-7 


peed 


TZ 


TPl 


TP  2 


coap 

calc 

coap 

calc 

coap 

calc 

coap 

0738 

22.57 

22.95 

12.946 

12.946 

107.583 

107.583 

1061 

15.67 

15.70 

9.014 

9.006 

75.  130 

74.846 

1477 

11.28 

11.28 

6.470 

6.467 

53.793 

53.793 

Analytical  equations  used  to  calculate  seccnd-crder 
Roaotc  transfer  function  coefficients  are: 

K  =H  /N  ;  T  =-  (I  -N .  )  /N 

5  r  2  r  r 

Analytical  equations  used  to  calculate  third-crder 

transfer  function  coefficients  are: 

K  =  <H6  *  VW /(W 

TZ  =-{  (B-V  *S6  -H^  )  /{IV*S6  '  V*5  > 

TP  1  ♦  IP  2=  (  (H- Y^)  *Hr  ♦  ( I2  -  )  *  Yv  ♦  *  (  Yr  - 8*  0)  ♦  Y^;  *  Nv  ) 

/(V( 

The  nondiaensionalized  hydrodynaaic  coefficients  for  the 
SL-7  containership  are  shown  in  Table  2. 


TABLE  2 

HYDBODYNAHIC  COEFFICIEHTS  FOB  THE  SL-7 


axial 

force 

lateral 

force 

aooent 

z-axis 

X  *  • 

u 

*  -0.0001 

l'v  * 

-0.00758 

H’v  " 

-0.00213 

X* 

uu 

=  -0.0003 

x,r  = 

0.0023 

H*  = 

r 

-0.00105 

X* 

vr 

*  0.0039 

0.00145 

H»  * 

6 

-0.0007 

*'vv 

=  -0.0012 

Y«  * 

vvr 

0.01 

N' 

vvr 

-0.015 

X*66 

*  -0.0005 

Y»  = 

vrr 

-0.008 

I* 

vrr 

-0.008 

Y' 

vw 

-0.03 

»* 

vw 

0.01 

N 

I*  s 

rrr 

0.003 

H*  » 

rrr 

-0.006 

▼  I  S 

M  »  (>  f 

-0.0005 

w  *  = 

0  »  »  • 

0.0001 

Ill,  £OST  FOHCTIOS 


10  recent  years,  oany  have  studied  the  problem  of 
£Bef.  2]  £Bef.  3]  £Bef.  4]  £Bef-  5]  £Bef.  *]  £Bef.  7] 
£Bef .  8]  £Bef .  9]  £ Bef.  10]  £Bef.  11]£Bef.  12]  £Bef.  13] 
£Bef.  14]  optimizing  an  automatic  ship-steering  controller 
for  minimum  fuel  consumption.  It  is  well  known  that  addi¬ 
tional  drag  is  introduced  by  steering  and  that  both  the 
rudder  motion  and  the  yawing  motion  contribute  to  this  added 
drag.  A  measure  of  the  added  drag  given  as  a  cost  function 
is 


i 

j  =  i/t^{  ♦  a  *  j 


dt 


(3-1) 


where  =  yaw  error 

5  =  rudder  angle 
X  =  weighting  factor 

tihile  this  expression  is  an  approximation,  it  is  conven¬ 
ient  for  shipboard  use  because  and  i  are  readily  measur¬ 
able.  There  is  no  general  agreement  on  numerical  values  for 
the  weighting  factor,  X  ,  and  in  this  study  the  values  used 
were  chosed  from  the  work  of  B.E.  Reid  [Be f.  7]  for  the 
SL-7. 

The  weighting  factors  for  the  operating  range  of  the 
ship  are  shown  in  Table  3. 


TABLE  3 

WEIGHTING  FACTOR 


ship  speed 
(knots) 


16 

23 

32 


weighting  factor 


16.796 

8.128 

4.2 


Reid's  work  shows  the  relationship  of  weighting  factor 
to  the  closed-loop  natural  frequency,  Bass,  pivot  point, 
ship  speed,  X"vr  and  X' hydrodynanic  coefficients.  It  is 
shown  in  Equation  3.2.  Reid  chose  a  closed-loop  natural 
frequency  of  0.05  rad/sec  which  experimentally  shewed  at 
this  frequency,  the  weighting  factor  in  the  cost  function, 
provided  good  representation  of  the  added  drag  due  to 
steering. 


2*fl*(1*X"  )*  (CP/L)*»*  /(P/2)  *(WX..  *0*) 


(3.2) 


XT.  . COlTiOLLM  DESIG8  OSIBg  ICSOS 

The  Interactive  Control  Systea  Optinization  and 
Sinulation  (ICSOS)  package  finds  optiaua  values  for  unknown 
(free)  paraaeters  in  a  control  sjstea  design  problea  and/or 
perfozas  siaulation  of  the  systea.  An  exaaple  of  nsage  of 
ICSOS  is  shorn  in  Appendix  B. 

In  preliainary  studies  ICSOS  was  used  with  Boaotc  nodels 
to  study  controller  characteristics  in  cala  water.  The  func¬ 
tion  ainiaization  subroutine  adjusted  the  controller  paraae¬ 
ters  to  liniaize  the  cost  function.  Figure  4.1  shows  the 
scheae  used  to  evaluate  the  controller  paraaeters. 


NOMOTO 

MOOEL 


Adjust 

Poromstsrs 


FUNCTION 

MINIMIZATION 

SUBROUTINE 

jW  +8*)dt 


Beid  [fief.  7]  uses  the  second~order  Nomoto  model  of 
equation  2.1  for  the  SL-7  and  also  uses  a  controller 
described  by 


Gc(S)  =  K 1* ( 1  *T1*S)  /  ( 1*T2*S) 


(4.1) 


Bis  results  are  given  in  Table  4. 


TABU  4 
fiEIO'S  BESULTS 


speed 

plant 

weighting 

controller 

gains 

knots 

K 

T 

factor 

K 1 

T 1 

T2 

16 

0.  1084 

90.36 

16.796 

0.4556 

89.51 

10.06 

23 

0.1556 

64.67 

8.128 

0.3769 

62.60 

8.308 

32 

0.2167 

45.45 

4.2 

0.3188 

44.92 

7.066 

Using  this  plant  and  weighting  factor  values  but  applying 
ICSOS ,  results  were  obtained  and  shown  on  Table  5. 


TABLE  5 
ICSOS  BESOLTS 


speed 

plant 

weighting  controller 

gains 

cost 

knots 

K 

T 

factor  K1 

T1 

T2 

J  min 

16 

.  1084 

90.36 

16.796  .454616 
8.128  .373171 

90.3459 

10.0215 

340.864 

23 

.1556 

64.67 

64.6658 

8.4640 

139.9916 

32 

.2167 

45.45 

4.2  .318645 

45.4475 

7.0662 

60.828 

In  each  case  the  controller  zero  (1/T1)  cancels  the 
plant  pole  (1/T).  Additional  experiments  consisting  of 
inserting  arbitrary  numbers  in  the  Boaoto  equation  and 
repeating  the  coaputei  run  indicated  that  this  will  always 
be  true.  That  is,  to  nininize  the  cost  the  plant  pole  is 
cancelled  and  a  new  pole  location  determined  with  appropri- 
ately  adjusted  gain. 


The  siaple  controller  of  Equation  4,1  is  an  arbitrarily 
chosen  structure.  To  deteraine  the  effects  of  aore  coaplex 
controllers  three  additional  structures  were  chosen  as  shown 
in  Figure  4.2.  Each  of  these  was  used  with  the  Hcaoto 
second- order  aodel  for  the  ship  at  each  of  the  indicated 
speeds.  The  results  are  shown  in  Tables  6,  7,  and  8. 


Figure  4.2  V1BI0US  STB  OCT  DEES  FOB  COBTBOUEBS 


These  results  are  very  interesting,  it  16  knots  the 
controller  gain  (K1)  ,  controller  zero  (1/T1)  and  ccntroller 
pole  (1/T2)  are  essentially  the  saae  for  all  structures.  For 
structure  B,  which  includes  an  additional  pole,  the  function 
ainiaization  subroutine  tries  to  drive  the  additional  pole 
to  infinity,  and  no  doubt  would  have  done  so  if  the  calcula¬ 
tions  had  continued.  For  structure  C,  which  has  two  poles 
and  two  zeros,  a  zero  and  pole  cancel  indicating  that  they 
are  not  needed  or  wanted.  For  structure  D,  the  integrator 


TABLE  6 

SIMULATION  RESULTS  -  STEADY  STATE  600  SECS 


CALB  BATES  FOB  VABIOOS  CONTROLLERS 
FOB  FIXED  SHIP  SPEED  (  16  KNOTS  ) 
10B0T0  SECOND-ORDER  BODEL  IK=. 1084 . T=90. 36) 
X  =  16.796  .  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 


contr 


A  . 

B  . 

C  . 

D  . 


K1 

454616 

<144101 

454511 

454581 


controller  gains 
T 1  T2  T3 

90.435S  10.0215 
90.2950  9.8566  .01 

90.3685  10.0224  23.  085 
90.37 1 S  10.0222 


T4 

23.084 


cost- 
Ti  J  ain 

340.864 
341.046 
340.864 
1E09  340.864 


TABLE  7 

SIMULATION  RESULTS  -  STEADY  STATE  600  SECS 


CALM  BATES  FOR  VARIOUS  CONTROLLERS 
FOB  FIXED  SHIP  SPEED  (  23  KNOTS  ) 
NOBOTO  SECOND-ORDER  MODEL  fK=. 1556. T=64. 67) 
X  =  8.128  .  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 


contr 


K1 


A  .373171 
B  .340024 
C  .373139 


controller  gains 
TI  T2 

64.6657S  8.463957 

79.65872  8.889204 

64.66855  8.463497  25 


T3  T4 

.01 

.9719  25.9738 


cost 
J  Bin 


139.9916 
140.9338 
13  9. 9991 


TABLE  8 

SIMULATION  RESULTS  -  STEADY  STATE  600  SECS 


CALM  BATES  FOB  VABIOOS  CONTROLLERS 
FOB  FIXED  SHIP  SPEED  f  32  KNOTS  1 
NOMOTO  SECOND-OBDER  MODEL  IK=. 21 67 .1=45. 45) 
X  =  4.2  .  OPTIMAL  PARAMETER  GAINS  OF 

VARIOUS  CONTROLLERS  ,  COST  FUNCTION 


contc  controller  gains 

K 1  TI  T2  T3  T4 

A  .318645  45.44747  7.06617 
B  .318  45.45  7.066  .05 

C  .318678  45.57511  7.06790  50.1829  50.04832 


cost 
J  ain 


60.828 

60.933 

60.828 


gain  is  driven  to  2ero.  The  saae  pattern  of  results  is 
obtained  at  23  knots  and  32  knots.  Note  that  in  all  cases 
the  liniaua  cost  is  essentially  the  sane,  as  vould  be 
expected  since  all  controllers  are  the  saae. 

Using  the  computer  aethod  of  Figure  4.1  and  the  Ncaoto 
third-order  node  Is  of  Table  1,  controllers  A,  B,  C  of  Figure 
4.2  vere  optimized..  The  results  are  shown  in  Tables  9,  10, 
and  1 1. 


TABLE  9 

SI  ISOLATION  BESDLTS  -  STEAOI  STATE 


600  SECS 


CALH  BA TER  FOR  VARIOUS  CONTROLLERS 
FOR  FIXED  SHIP  SPEED  (  16  KNOTS  ) 

HOHOTO  THIRD-ORDER  HODEL 
(K=. 073812,12=22. 5673. TP  1=12. 9458, TP2*107. 5853) 
X  -  16.796  ,  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 


contr 


K 1 


controller  gains  cost 

T 1  T2  T3  T4  J  Bin 


A  0.6446104 

l 


90.CS94 

_  84.  826 

1139  107.5782 


15.2771 2  -  -  370.4023 
15.78691  .24598  -  374.3808 
8.73520  12.9368  24.9676  369.9297 


TABLE  10 

SIMULATION  RESULTS  -  STEADT  STATE  600  SECS 


contr 


CALH  BATEB  FOR  VARIOUS  CONTROLLERS 
FOR  FIXED  SHIP  SPEED  (  23  KNOTS  ) 
lOMCTO  THIRD-ORDER  MODEL 
(K*. 1067.TZ=15. 675 .TP 1*9. 014, IP2-75. 13) 
X*  8.128  .  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 


K1 


controller  gains 
T 1  T2  .13 


T4 


cost 
J  ain 


A  0.5224258  63.13609  12.72212 
n  «u_q nnq 


_  -  -  152.2920 
0.5216467  64.93709  12.63218  .0505174  -  152.5333 
€.5001907  75.14852  6.527490  9.039928  18.260  152.2800 


Of  aajcr  interest  is  the  fact  that  the  difference  in 
"cost"  between  A,  B,  C  is  less  than  one  per  cent.  At  each 


TABLE  11 

SIHULATIOH  RESULTS  -  STEADY  STATE  600  SECS 


CALH  BATES  FOB  YABIOUS  COHTBOLLEBS 
FOB  FIX EE  SHIP  SPEED  l  32  KHOTS  ) 

ICHCTO  THIBD-OBDEB  HO DEL 
(I=.  14771- TZ  =  1 1.2833.  TP1=6.4699,  TP  2=53.  7931) 
A  =  4.2  *  OPTIHAL  PABAHETEfi  GAIHS  OF 

YABIOUS  COHTBOLLEBS  ,  COST  FUHCTIOH 


contr 


A 

B 

C 


K1 

0.427633 

0.298732 

0.417991 


controller 
T 1 


chains 


48.66C48 

89.40696 

53.69961 


10.744  65 
15.01033 
4.970016 


T3 


.0597786 

6.294354 


T4 


13.85724 


cost 
J  Bin 

68.09039 
69.323  55 
68.04735 


speed  (16,23/32  knots)  controller  C  is  "BEST"/  tut  the 
difference  is  slight.  Examining  the  parameter  values 
obtained  fcr  controller  C,  it  is  seen  that  at  all  three 
speeds  both  poles  of  the  ship  are  essentially  cancelled  by 
zeros  of  the  controller. 

These  results  seem  to  indicate  that  the  dynaaics  cf  the 
plant  determines  the  cptiaua  structure  for  the  controller. 

Using  a  state-feedback  controller  and  Noaoto  third-order 
nodels  of  Table  1,  the  controller  was  optimized  for  various 
ship  speeds.  Figure  4.3  shows  the  schene  used  to  evaluate 
the  state-feedback  controller. 

Using  the  scheme  of  Figure  2.2,  with  no  change  in  cost 
function  or  weighting,  the  optimal  gains  and  costs  were 
deter lined  as  shown  in  Table  12. 

Hhen  comparing  the  state-feedback  controller  with 
controller  C,  it  is  seen  that  at  each  speed  controller  C  has 
a  lower  cost.  Among  the  controllers  consired,  controller  C 
is  "BEST"  when  using  the  Nonoto  third-order  model,  although 
the  differences  in  cost  are  not  dramatic. 


Tig are  5.1  0PTIBIZ11Z0H  OF  C0HTB01LBH  OSIMG  FOBTBAI  PBOGBAfl 


Using  the  coaputec  method  of  Figure  5. 1  and  the  nonli¬ 
near  equations  of  notion,  controllers  A,  B,  C  of  Figure  4.2 
vere  optimized.  The  results  are  shown  in  Tables  13,  14,  and 
15. 


TABLE  13 

SIBOLATIOH  BESOLTS  -  STEAD!  STATE  600  SECS 


CALB  BATEB  FOB  TABI00S  COBTBOLLEBS 
FOB  FIXES  SHIP  SPEED  (  16  KNOTS  ) 
EQUATIONS  OF  NOTION 

A  =  16.796  .  OPTIMAL  PABABETEB  GAINS  OF 
YABIODS  CCNTBOLLEBS  ,  COST  FONCTION 

contr  controller  gains 

K1  T 1  T2  T3  T4 

A  .648401  89.81704  15.381699 

B  .620050  90.67294  15.542297  0.9201336 

C  .617326  107.  1494  8.597198  13.353928  25.21362 


cost 
J  nin 


1.  128189 
1.  173323 
1. 126307 


TABLE  14 

SIBOLATIOH  BESOLTS  -  STEAD!  STATE  600  SECS 


contr 


A 

B 

C 


CALB  BATES  FOB  VABIOOS  CONTBOLLEBS 
FOB  PIXEL  SHIP  SPEED  (  23  KNOTS  ) 
EQUATIONS  OF  NOTION 

A-  8.128  .  OPTIMAL  PABABETEB  GAINS  OF 
VABIOOS  CCNTBOLLEBS  ,  COST  FUNCTION 


K1 


T 1 


controller  gains 
T2  T3 


522106  66.33122 
,495869  66.15152 
503967  74.79771 


12.83327 

13.01183 

6.65880 


0.92783 

9.20533 


T4 


18.4022064 


cost 
J  min 

0.4640879 

0.4857854 

0.4636095 


These  results  agree  with  those  obtained  by  ICSOS  and 
controller  C  provides  the  minimum  cost. 

If  the  assumption  that  the  steering  dynamics  of  the  ship 
is  adequately  modeled  as  a  second-order  system  is  valid, 
then  only  two  states  are  needed  for  feedback.  For  a  third- 
order  system  three  states  are  required.  The  controller 
structures  are  shown  on  Figure  5.2  and  5.3.  Using  the  scheme 
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TABLE  15 

SIHOLATIOB  RESULTS  -  STEAD!  STATE  600  SECS 


m 


tea 


t  - 


contr 


K1 


CALS  RATER  FOR  VARIOUS  COBTROLLERS 
FOB  FIXED  SHIP  SPEED  (  32  KBOTS  ) 

EQUATIOBS  OF  MOTIOB 
»  4.2  .  OPTIMAL  PARAMETER  GAIBS  OF 

VARIOUS  COBTROLLERS  ,  COST  FUBCTIOB 

controller  gains  cost 

T1  T2  T3  T4  J  nin 


A 

B 

C 


.428404  48.65540  10.814426 
.298732  89.40696  15.010330 
.417333  “  -  - 


53.0965  4 


0.2072417 

_  0.01  -  0.2118334 

5.096548  6.474857  14.0205  0.2071124 


of  Figure  5.1  with  no  change  in  cost  function  or  weighting, 
the  optiaal  gains  and  cost  uere  determined  as  shown  in 
Tables  16  and  17.  Coiparing  costs,  there  is  little  differ¬ 
ence  between  the  twc  state  systea  and  the  three  state 
systea.  The  coaparison  between  state  feedback  with 
controller  C,  it  is  seen  that  at  each  speed  controller  C  is 
better,  but  not  auch  better. 


Figure  S.3  THEBE  STUB  SISTEH 


TABLE  16 

I 

SIHOLATIOH  HISOLTS  -  STEADI  STATE  600  SECS 


speed 

knots 


/ 


CALE  VATIB  FOB  YABIOOS  SHIP  SPEEDS 
2Q0ATI0BS  OF  HOTIOH 
OPTIBAL  PABABETEB  GAIHS  FOB 
1BO  STATE  SISTEH 


K1 


gains 


K2 


4.4033689 

3.0889006 

2.2342062 


mm 


weighting  cost 
factor  J  ain 


16.796  1.128771 

8.128  .4646050 

4.2  .2075207 


Bote  that  for  the  Bonoto  model  studies  yaw  error  and 
rudder  angles  were  Measured  in  degrees;  when  the  eguations 
of  notion  were  siaulated  yaw  error  and  rudder  were  in 
radians.  Thus  the  numerical  values  of  the  cost,  J,  are 
different.  y 

Transient  response  plots  for  controllers  A,  B,  C,  and 
three  state- feedback  at  ship  speed  32  knots  are  shown ‘in 
figures  5.4  -  5.9. 
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TABLE  17 

SIBOLATIOB  BZSOLTS  -  STEAOZ  STATE  600  SECS 


speed 

knots 


4. 


CALB  BATES  FOB  VABXOOS  SHIP  SPEEDS 
EQOATIOBS  OP  HOTIOB 
OPTIBlL  PABABETEB  6AIHS  FOB 
TBBEE  STATE  SISTEH 


K1 

>617249 


gains 

*K2 

K3 

'liltH19 

cost 

J  *in 

87.7073364 

56.2784882 

33.7511444 

99.9802704 

88.5913391 

41.3186035 

16.796 

8.128 

4.2 

1.  128289 
.464  3548 
.2074225 

32  KNOTS 


TRW  VS  TlME-FEEDSflCK  flNO  fl.B.C  COMP 


32  KNOTS 

YftHUOOCR  VS  TlfC-COMPENSfllOR  B 


T«N 


RUOOCR 


71.  C0H1E0LLEB  DESIGN  IN  SEA  STATE 


A  sea  state  generator  was  coupled  to  the  Fortran 
program,  so  that  the  function  minimization  subroutine  could 
be  used  in  the  presence  of  the  sea  state.  The  sea  state 
generator  was  an  elaborate  program  obtained  from  DTNSBDC. 
This  program  generates  added  mass  and  inertia  values  as 
functions  of  encounter  frequency  and  also  calculates  the 
forces  and  moments.  The  forces  and  moments  are  generated 
and  stored  in  a  look  up  table  which  was  coupled  to  the  equa¬ 
tions  of  motion.  Figure  6. 1  shows  the  scheme  used  to  eval¬ 
uate  the  controller  parameters.  The  computer  program  is 
shown  in  Appendix  A. 

The  optimal  gains  obtained  by  the  calm  water  study  of 
Chapter  five  were  use  as  the  initial  guess  in  evaluating  the 
optimal  controller  parameters  in  the  presence  of  a  seaway. 
For  comparison,  studies  of  the  value  of  the  cost  function 
using  calm  water  gains  in  sea  state  were  obtained;  then  the 
function  minimization  subroutine  was  allowed  to  adjust 
controller  parameters  in  the  presence  of  several  sea  states 
and  encounter  angles.  The  entire  study  was  done  at  a  ship 
speed  of  32  knots.  The  added  mass  and  inertia  change  with 
respect  to  encounter  frequency  as  shown  in  Figures  6.2  and 
6.3.  Figure  6.3  is  ncndimensionalized  by  dividing  the  added 
inertia  by  the  mass  of  the  displaced  water  and  the  square  of 
the  length  between  perpendiculars.  To  convert  back  to  dimen- 
sionalized  units  of  lb-ft-sec2,  multiply  the  graph  points  by 
2.581112.  Since  the  sea  state  is  represented  by  irregular 
waves,  the  waves  impinging  on  the  ship  hull  contain  the 
total  energy  density  spectrum  composed  of  many  frequencies 
and  the  ship  responds  to  an  average  value  of  added  mass  and 
inertia.  The  values  used  for  this  study  was  obtained  at 


encounter  frequency  cf  0.75  rad/sec  from  our  sea  state 
generator.  This  frequency  gave  us  values  for  added  Bass  and 
inertia  representative  of  an  average  value.  The  energy 
density  spectra  for  various  sea  states  are  shown  in  Figure 
6.4.  The  added  nass  for  sway  was  changed  from  2.6457E06 
lb-sec2/ft  for  cala  water  to  2.3043EQ6  lb-sec2/ft  for  a 
seaway.  The  added  aoaent  of  inertia  for  yaw  was  changed 
froa  1.42E11  lb-ft-sec2  for  cala  water  to  1.5096111 
lb-ft-sec2  for  a  seaway.  All  other  hydrodynamic  coeffi¬ 
cients  were  kept  constant  at  cala  water  values.  The  results 
are  shown  in  Tables  18  -  25.  In  certain  sea  states  and 
encounter  angles  the  cala  water  optiaal  gains  performed  well 
as  shewn  by  cala  water  cost  value  when  conpared  to  sea  state 
cost  value.  In  most  cases  the  function  miniaization  subrou¬ 
tine  found  new  gains  with  lower  cost  function  values  in 
seaway  as  compared  to  using  cala  water  gains.  In  the  cala 
water  evaluation,  the  system  was  perturbed  with  a  one  degree 
course  change,  but  tke  course  change  was  not  included  in  the 
seaway  tests.  The  difference  in  cost  values  is  attributed  to 
the  difference  in  operating  conditions. 

Osing  the  Proportional,  Integral  and  Derivative  (PID) 
controller  Eguation  €.1  with  no  change  in  cost  function  , 
the  function  minimization  subroutine  was  used  to  adjust 
controller  parameters  to  miniaize  the  cost  function  and 
evaluate  the  minimum  cost.  The  results  are  shown  in  Table 
26.  Bhen  comparing  the  PID  with  controller  A,  it  is  seen 
that  at  each  encounter  angle,  controller  A  is  better.  These 
results  agree  that  in  a  seaway  controller  A  provides  the 
minimum  cost. 

Table  27  shows  comparison  of  the  minimum  cost  function 
for  controller  A,  controller  C,  and  PID.  The  study  was  done 
at  ship  speed  of  32  knots  and  at  sea  state  4.  Controller  A 
provides  the  minimum  cost. 


Figure  6.1  OPTIHIZATIOH  OF  COHTBOL1EB  III  SEA  STATE 


The  optimal  gains  obtained  in  the  presence  of  sea  state 
was  done  over  using  a  simulating  tine  of  600  seconds.  The 
sea  state  program  is  designed  to  provide  gradual  increase  in 
the  forces  and  moments  during  an  initial  time  interval. 
This  is  done  to  minimize  initial  condition  transients  in  the 
ship  dynamics.  There  willr-  unavoidably  be  some  transient 
effects,  however,  and  these  could  affect  the  value  of  the 
cost,  J,  determined  during  the  600  seconds  of  simulation.  To 
determine  whether  such  initial  transients  had-  any  signifi- 
cant  contribution  to  the  value  of  J,  additional  simulation 
runs  were  made  with  the  controller  parameters  fired  at  their 
optimal  values.  However,  evaluation  of  the  cost,  J,  was 
started  only  after  300  seconds  of  simulation  had  elapsed. 


ADDED  MASS  (*10"6) 


ADDED  MASS  COEFFICIENTS 

SWAY  WRT  SWAY  ACCELERATION 


Figure  6.2  AEDZD  BASS  vs.  BHCOONTEH  FREQOEHCI 
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NONDIMENSIONAL  COEFFICIENTS 

0.00  0.07  0.14  0.21  0.28  0.35  0.42 


ENCOUNTER  FREQUENCY  (RAD/SEC) 


figure  6.3  AOSZD  IHEBTIA  vs.  EHCOONTEB  FBBQOEHCI 
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AMPLITUDE  SPECTRA  (FT"2-SEC) 


TABLE  18 

SIHOLATION  BESOLTS  -  STEADY  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KHOTS)  IN  A  SEABAY 
SHIP  HODEL:  EQUATIONS  OF  HOTIOI 
YAHC-0.0 

SEA  STATE  1 


COHTBOLLEB  A 


encounter  controller  gains 

sea  state 

cost  with 

angle 

cost 

calm  water 

degree 

K1 

11 

72 

J  ain 

J 

0 

.4284037 

48.6554395 

10.814426 

. 6 17452- 34 

. 61745E-34 

30 

1. 1561117 

29.3693695 

1.4592390 

.2870198 

.5128402 

60 

1.4033298 

10.6530075 

1.  1086683 

. 1342071 

.2154726 

90 

.2969198 

58.2413940 

1.8758221 

. 1300669 

.1565958 

120 

.  1761794 

299.999512 

30.7967834 

.05741726 

.0727727 

150 

2.8430557 

5.2826872 

.8887696 

.0219070 

.0939  400 

180 

1.6211386 

14.0782928 

2.0712433 

.0051925 

.0095694 

TABLE  19 

SIHOLATIOH  BESOLTS  -  STEADY  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEABAY 
SHIP  HODEL:  EQUATIONS  OF  HOTION 
YABC-0.0 

SEA  STATE  2 

CONTBOLLEB  A 


encounter 

angle 

degree 


controller  gains 


K1 


IS 

90 

120 

150 

180 


.42840370 
.27997030 
.95575100 
1.3577642 
1.  1208973 
2.S777727 
.61420630 


II 

48.6554395 
249.935059 
24. 3813629 
9.49564080 
25.4498596 
16.2154541 
.482041200 


T2 

10.814426 

1.  1068363 
4.0224676 
.56274800 
6.2521963 


sea  state 
cost 
J  ain 

.617452-34 

.04774852 

.04104504 

.02650556 

.04928402 

7.5751530 

.000124338 


cost  with 
calm  water 
J 

. 61745E-34 

. 0886  22  5 

.0535875 

.0483197 

.0717524 

28.1294403 

.0002445 


TABLE  20 

SIBOLATIOH  BESOLTS  -  STEADY  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KHOTS)  IN  A  SEA8AY 
SHIP  HODEL:  EQOATIOHS  OF  HOTIOH 
YA1C=0. 0 

SEA  STATE  4 


COHTBOLLEH  A 


encounter  controller  gains 

sea  state 

cost  with 

angle 

degree 

K1 

T1 

T2 

cost 

J  Bin 

cals  water 

J 

,0  ’ 

.4284037 

48.£££40 

5.7*3536 

10.814426 

. 620598E-34 
.02854677 

. 62C598E-34 

30 

.9815440 

. 6999879 

.0395892 

60 

.6201209 

40.  8C556 

19.606873 

.09375697 

-1032696 

90 

1.809746 

36. 01225 

6.324708 

1.5171340 

4.  1623011 

120 

5. 195190 

18. 92513 

.  6999907 

9.991730 

48.970703 

150 

1.446776 

16.89375 

. 5265408 

16.67052 

24.822098 

180 

.1000000 

1.000000 

20. 149999 

.00739631 

-0076657 

TABLE  21 

SIBOLATIOH  BESOLTS  -  STEADY  STATE  600  SECS 

EIXED  SHIP  SPEED  (32  KHOTS)  IH  A  SEAHAY 
SHIP  HODEL:  EQUATIOHS  OF  BOTIOH 
YAVCs0. 0 

SEA  STATE  6 

COHTBOLLEB  A 


encounter 

angle 

degree  K1 


controller  gains 

T 1  T2 


o 

30  2.5715786 

60  1.7228041 

90  1.6584366 

120  3.3422489 

150  .2854474 

180  . eC53379 


48.6  ££3955 

10.4721832 

8.4014740 

37.1672655 

106.722259 

157.483887 

.75733550 


10.8144264 
. 5342450 
. 5141125 
.5792384 
. 9260592 
119.981018 
6.04484460 


sea  state 
cost 
J  Bin 

50899E-32 
1.4287940 
1.5827220 
4.5505371 
22. 108002 
.81100580 
.07365978 


cost  with 
call  water 
J 

.  50899e- 32 
4.747240  10 
3.427449  20 
13.2757149 
94.5497589 
1.50448510 
.142564400 


TABU  22 

SIHOLATIOB  BESOLTS  -  STEADI  STATE  600  SECS 


EIZED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAHA1 
SHIP  HCDEL:  EQUATIONS  OF  HOTION 
I INC =0.0 

SEA  STATE  1 

CONTBOILEB  C 


encounter 
angle 
degree  K1 


controller  gains 
T1  T2  T3 


30 

60 

90 

120 

150 

180 


;  5.0481 


1.61345  16.  8755 
.957558  12.6178  7.32113 
.781984  17.6475  9.22485 
.417332  53.0965  5.09655 
.417332  53.0965  5.09655 
2.  13735  18.8265  17.5778 
.957558  12.6178  7.32113 


47.3405 
43.3531 
13.9438 
6.47857 
6.47857 
25. 1516 
43.3531 


T4 

1.73759 
8.15752 
16.76  63 
14.0205 
14.0205 
21.1481 
8.  15752 


sea 
cost 
J  oin 

,  14E-33 
.324739 
.  159710 
,  148588 
,077918 
,031496 
.006566 


ca  1b 
cost 
J 

.4 
.3 
.203137 
.148588 
.077918 
.081612 
.008  172 


TABLE  23 

SIHOLATIOB  BESOLTS  -  STEADI  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAHAI 
SHIP  HODEL:  EQUATIONS  OF  HOTION 
TAHCsO.O 

SEA  STATE  2 

CONTBOILEB  C 


encounter 
angle 
degree  K1 


controller  gains 
T1  T2  T3 


0 

to 

90 

120 


18 


.417333 

.849594 

.417333 

.781984 

.880395 

.899999 

.440916 


sea 
cost 
T4  J  Bin 


5  3  .  _  _  . 
i9.99i; 
53.0961 
17.6475 
21.4597 
15.7103 
.093671 


5.09655 

9.34138 

5.09655 

9.22485 

10.9255 

1.11632 

17.7305 


6.47486 


20.3578 

.47486 


13.9438 

9.24547 

41.3275 

25.2103 


14.0205 

13.7487 

14.0205 

16.9438 

11.0667 

2.29012 

5.04178 


caln 

cost 

J 


. 18E-33 
.054338 
.044536 
.033518 
.055292 


10,763 

.00012 


.  18E-33 
-061184 
.044536 
.048467 
.0  56288 
23.8522 
.001635 


TABLE  26 

SIHOLATION  BESOLTS  -  STEADY  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KHOTS)  IN  A  SEAWAY 
SHIP  HCDEL:  EQUATIONS  OF  NOTION 
T ABC -0.0 


SEA  STATE  4 
PID  CONTBOLLEB 


encounter 

controller  gains 

sea  state 

angle 

K1 

T1 

T2 

J  ain 

30 

.95263100 

4.20720860 

.69368610 

.029  85619 

60 

.68631890 

12.5794449 

8.2121658 

.09730512 

90 

2.5809155 

12.4247589 

.77810380 

1.5915950 

120 

4.9198265 

12.5986176 

.67592390 

10.  708980 

150 

1-3970823 

15.7682953 

.51991  180 

17.427200 

8(S)/^(S)  = 


K 1  ♦  K1*T1*S  /  (  1  *T2*S)  **2 


(6.1) 


TABLE  27 

COSPABISON  OF  TEE  HIHIflUH  COST 


encc outer 
angle 
degree 

30 

60 

90 

120 

150 


SHIP  SPEED  (32  KHOTS) 
SEA  STATE  4 


controller 

A 

J  ain 

.02854677 

.09375657 

1.5171340 

9.9917300 

16.670520 


controller 

C 

J  ain 

.034033 
.098914 
1. 57368 
10.3530 
20.3956 


controller 
PID 
J  ain 

.02985619 

.09730512 

1.5915950 

10.708980 

17.427200 


The  value  obtained  was  then  doubled  and  conpared  with  the 
result  of  evaluating  3  over  the  full  600  seconds.  Coaparison 
of  Table  28  with  cost  values  in  Tables  18,  19,  20,  21  shows 
only  snail  differences. 

To  obtain  insight  into  the  stochastic  process  of  irreg¬ 
ular  seas,  a  deterninistic  process  was  studied.  The  Fortran 


TABLE  28 

EFFECTS  DOE  TO  THANSIENT  ADD  GHADOAL  BOILD  OP  OF  SEA  STATE 


IITEGBATIOI  OF  COST  FUNCTION  (  300  TO  600  SECS) 
FIXED  SHIP  SPEED  (32  KHOTS)  IN  A  SEABAI 


sea  encounter 
state  angle  K1 


controller  gains  cost  cost 

T 1  T2  J  Bin  2*J  Bin 


1  60  1 .403329 fi  10.650075  1.  1086683  .0641  122  .12-82244 

2  60  . 95575  ICO  24.381363  2.  3079853  .0199731  .0399462 
4  60  .62012090  40.805560  19.606873  .0515974  .1031948 
6  60  1.7228041  8.4014740  .51411250  .7906179  1.58123b 


program  was  modified  to  minimize  the  cost  function  in  the 
presence  of  a  regular  sea.  To  allow  comparison  with 
previous  work  the  encounter  frequency  of  0.75  rad/sec  was 
used  and  scaled  the  amplitude  of  the  regular  sea  to  its 
prospective  sea  state.  The  entire  study  was  done  at  a  ship 
speed  of  32  knots.  The  results  are  shown  in  Tables  29 ,  30, 
and  3  1. 

Table  29  shows  that  for  regular  seas  the  controller 
parameters  do  not  change  significantly  for  different  sea 
states;  but  as  sea  state  increases,  the  cost  value  increases 
due  to  the  increase  in  yaw  aoaent  and  sway  force  on  the 
ship.  Tables  30  and  31  also  show  that  the  controller  parame¬ 
ters  do  not  change  significantly  froa  sea  state  to  sea 
state.  However,  an  encounter  angle  of  90  degrees  shows  a 
relatively  high  cost  compared  to  costs  calculated  for  60  and 
120  degrees  at  a  given  sea  state.  To  account  for  this 
anomaly,  the  following  is  suggested.  In  the  regular  sea, 
the  added  aass  and  inertia  were  known  for  a  given  encounter 
frequency,  while  in  the  irregular  sea  a  representive  average 
value  was  used.  The  method  used  to  obtain  the  average  might 
not  represent  the  actual  average.  Also,  it  seeas  reasonable 
to  suppose,  that  the  assumptions  of  the  function  weighting 
factor  are  satisfied  for  all  encounter  angles;  that  is,  the 
weighting  function  (Eg.  3.2),  which  appears  in  the  cost 


function  (Eg.  3.1}#  does  totally  represent  the  added  drag 
for  all  encounter  angles.  Future  study  is  needed  tc  answer 
these  guestions. 

The  sea  state  in  the  deterainistic  aodel  is  represented 
by  regular  waves.  On  this  description#  the  waves  iapinging 
on  the  ship  hull  correspond  to  only  one  freguency  in  the 
energy  density  spectrua.  In  the  case  of  irregular  seas# 
however#  the  spectral  coaponents  change  for  different 
states#  as  shown  in  Figure  6.4.  Thus  coaparison  of  the 
controller  paraaeters  obtained  for  regular  seas  with  results 
for  irregular  seas  is  not  justified.  The  function  ninisiza- 
tion  subroutine  adjusted  controller  paraaeters  to  niniaize 
the  cost  function  for  either  case  (irregular  or  regular 
seas)  as  shown  in  tables  32  and  33. 

TABLE  29 

SIBOLATIOV  B1S0LTS  -  STEADY  STATE  600  SECS 


FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  BEGOLAB  SEA BAY 
SHIP  BCEEL:  EQUATIONS  OF  NOTION 
YA1C=0.0 

ENC00NTEB  FBBQOBNCX  =  0.75  BAD/SEC 
EBCOUBTEE  ANGLE  =60  DEGBEES 


sea 

state 


COITBOLLEB  A 


controller  gains 
T1 


1449795 

1534657 

1514665 

1533340 


141.383179 

129.987473 

135.798737 

135.488495 


32.9405670 
3  1.4042358 
3 2.9749756 
33. 5585632 


cost 
J  Bin 

,000764582 

.003056434 

.009345479 

.022174600 


Bote  that  in  bcth  the  deterainistic  and  stochastic 
aodels#  aaong  the  controllers  considered#  controller  A  is 
"BEST"  in  a  seaway  disturbance#  although  the  differences  in 
cost  are  not  draaatic. 

Finally#  the  observed  dependence  of  optiaal  controller 
gains  on  sea  state  and  encounter  angle  suggests  that  an 


TABLE  30 

SIMULATION  RESULTS  -  STEAD?  STATE  600  SECS 


FI1ED  SHIP  SPEED  (32  KNOTS)  IN  A  REGULAR  SEABA? 
SHIP  HCDEL:  EQUATIONS  OF  MOTION 
TABC-0.0 

ENCOUNTER  FREQUENCY  =  0.75  BAD/SEC 
SEA  STATE  4 


CONTROLLER  A 


encounter 

controller  gains 

cost 

angle 

K  1 

T 1 

T2 

J  min 

30 

.2351043 

102.021973 

28.3396912 

.002985300 

60 

.  1514665 

135.7987  37 

32.9749756 

.009345479 

90 

.4964442 

66.546493 

49.7598267 

.048143090 

120 

.  1327230 

149.540543 

33.6013489 

.038937880 

150 

.4536914 

70.566528 

31.5839539 

.062534153 

TABLE  31 

SIMULATION 

RESULTS  -  STEADY 

STATE  600 

SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  REGULAR  SEAHAI 
SHIP  HCDEL:  EQUATIONS  OF  MOTION 
YABC=0. 0 

ENCOUNTER  FREQUENCY  =  0.75  RAD/SEC 
SEA  STATE  6 


CONTROLLER  A 


encounter 

controller  gains 

cost 

angle 

K 1 

T 1 

T2 

J  min 

30 

.2370022 

100. 122940 

28.0581207 

.007092878 

60 

.  1533340 

135.488495 

33.5585632 

.022174600 

90 

.5210407 

62.153702 

49.9858093 

.112772880 

120 

.  1414837 

142.695160 

35.3171234 

.091541650 

150 

.4587426 

71.451385 

33.4568024 

.  144615829 

adaptive 

controller 

■ust  be  used 

to  provide 

a  continuous 

■ini bub  on  the  cost  function. 


After  obtaining  the  optimal  gains  for  controller  A ,  to 
observe  the  behavior  of  the  rudder  and  yaw  motion  of  the 
ship,  transient  response  plots  were  obtained  for  controller 
A  at  ship  speed  of  32  knots  and  sea  state  4  for  various 
encounter  angles  as  shown  in  Figures  6.5  -  6.14.  Note  the 


TABLE  32 

COHPABISCN  OF  IBBEGDLAB  TO  BEGULAB  SEAS  COHTBOLLEB  GAIHS 


SEA  STATE 

4 

COHTBOLLEB 

A 

ncounter 

controller  gains 

angle 

K 1 

T1 

T2 

30 

(irregular) 

.9815440 

5.733036 

.6999879 

30 

(regular) 

.2351043 

102.021973 

28.3396912 

60 

(irregular) 

.6201209 

40.805560 

19.6068730 

60  • 

(regular) 

.  1514665 

135. 798737 

32.9749756 

90 

(irregular) 

1.809746 

36.012250 

6.3247080 

90 

(regular) 

. 4964442 

66.546493 

49.7598267 

120 

(irregular) 

5.195190 

18.925130 

.6999907 

120 

(regular) 

.  1327230 

149.540543 

33.60134  89 

150 

(irregular) 

1.446776 

16.893750 

. 52654C3C0 

150 

(regular) 

.4536914 

70.566528 

31. 5839539 

TABLE  33 

COHPABISCN  OF  IBBEGOLAB  TO  BEGOLAB  SEAS  COHTBOLLEB  GAINS 

SEA  STATE  6 


COHTBOLLEB  A 


ncounter 

controller  gains 

angle 

K  1 

T1 

T2 

30 

(irregular) 

2.9715786 

10-4721832 

.534  2450 

30 

(regular) 

.23700220 

100. 122940 

28.058  12C7 

60 

(irregular) 

1.7228041 

8.4014740 

.5141125 

60 

(regular) 

.  15333400 

135.488495 

33.5585632 

90 

(irregular) 

1.8584366 

37.1672655 

.5792384 

90 

(regular) 

. 5210407 

62.153702 

49.9858093 

120 

(irregular) 

3.3422489 

106.722259 

.9260592 

120 

(regular) 

.  1414837 

142.695160 

35.3171234 

150 

(irregular) 

.2854474 

157.483887 

119.981018 

150 

(regular) 

.4587426 

71.451385 

3  3.4568024 

increase  in  both  rudder  and  yaw  amplitude  as  the  encounter 
angle  increased.  This  is  due  to  the  increase  in  yav  moment 
and  svay  force  on  the  ship. 


32  KNOTS- SEfl  STATE  4 


ENCOUNTER  ANGLE  30  DEGREES 


32  KNOTS-SEA  STATE  i 
ENCOUNTER  ANGLE  30  DEGREES 


Figure  6.8  BDDDEB  vs.  TIME  60  0E6BEE5 
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32  KNOTS-SEfl  STATE  4 


ENCOUNTER  ANGLE  120  DEGREES 


i  T.  ■  i  i  .  . ■(  t  — —  . . . .  . . .  1  i" . i"  ■■  1  'i 

50.0  100.0  1S0.C  2CO.O  250.0  300.0  350.0  100.0  150.0  500. C  550.0  60 

TINE  " 


Figure  6.12  BOODEB  vs.  TIBS  120  DEGBEES 


VII-  IN  ADAPTIVE  CONTROLLER 


In  a  seaway,  the  controller  gains  changed  dramatically 
for  changes  in  sea  state  and  encounter  angle.  An  adaptive 
controller  oust  be  used  to  provide  continuous  operation  on  a 
minimum  of  the  cost  function.  This  Chapter  addresses  a 
theoretical  design  of  an  adaptive  controller. 

In  the  future,  there  will  be  better  measurement  of  navi¬ 
gation  than  can  be  provided  by  conventional  equipment  on 
board  a  ship.  Presently  the  Navy  is  involved  in  a  program 
that  will  provide  precision  navigation  data.  The 
NAVSTAH/GICBAL  POSITION  SYSTEM  (GPS)  [fief.  15]  [Hef.  16] 
[fief.  17]  will  provide  extremely  accurate  three-dimensional 
position  and  velocity  information  to  users  anywhere  in  the 
world.  The  position  determinations  are  based  on  the  measure¬ 
ment  of  the  transit  time  of  BF  signals  from  four  satellites 
of  a  total  constellation  of  eighteen.  This  system  is  sched¬ 
uled  to  be  fully  operational  in  1988.  At  present  (1984) 
there  are  four  N AVSTAB/GPS  satellites  in  operation  which 
allows  three  to  four  hours  per  day  of  navigation  time. 
Already  the  Texas  Instrument  Company  markets  a  receiver  for 
this  system  where  GPS  can  be  used. 

The  Navy  Remote  Ocean  Sensing  System  (N50SS)  [fief.  18] 
will  be  able  to  determine  wind  velocities  over  the  world's 
oceans  with  an  accuracy  sufficient  to  determine  ccean 
surface  waves.  It's  objective  will  be  to  acguire  global 
ocean  data  for  operation  and  research  use  by  both  the  mili¬ 
tary  and  civil  sectors.  This  system  is  scheduled  to  launch 
its  first  satellite  in  June  1989. 

The  scheme  for  an  adaptive  controller  is  shown  in  Figure 
7.1.  Having  stored  the  optimal  controller  parameters  in  a 
look  up  table  as  functions  of  ship  speed,  sea  state,  and 


encounter  angle,  the  ship  operating  condition  must  he  known 
so  that  the  table  is  useful.  KAVSTAB/GPS  would  identify 
ship  speed  and  NBOSS  would  identify  sea  state  and  encounter 
angle.  The  optimal  parameters  can  then  be  looked  up  and 
inserted  into  the  controller.  This  should  place  system  oper¬ 
ation  near  the  minimus  J.  To  ensure  fine  tuning,  a  micro¬ 
processor  programmed,  with  the  function  minimization  or-line 
in  machine  language,  with  inputs  of  yaw  error  and  rudder 
moticn  of  the  ship  would  accomplish  the  fine  tuning  rapidly. 
Since  the  subroutine  is  written  in  Fortran  (as  used  for  this 
study)  this  would  be  inappropriate  for  on-line  use. 

The  adaptive  controller  can  be  performed  with  digital 
circuits  rather  than  analog  components.  Garcia  [Bef.  19] 
demonstrates  the  process  for  converting  an  analog  controller 
into  a  digital  controller.  Figure  7.2  illustrates  the 
processing  of  the  major  components  in  a  digital  controller. 
An  analog  component  circuit  can  be  replaced  by  an  analog  to 
digital  converter,  a  digital  processor,  and  a  digital  to 
analog  converter.  Some  of  the  benefits  which  can  be  realized 
by  doing  this  are: 

1.  A  high-speed  processor  could  actually  process  a 
number  of  multiplexed  signals,  performing  processing  func¬ 
tions  on  a  number  of  independent  channels. 

2.  The  processing  function  is  permanent  in  software, 
unless  deliberately  changed,  and  will  not  drift  with  age. 

3.  The  processing  function  can  be  changed  without 
changing  components,  merely  by  changing  software. 

4.  Accuracy  can  be  made  very  high  and  can  be  changed 
merely  by  changing  software. 

5.  Processing,  which  previously  reguired  large  compo¬ 
nents  such  as  inductors  in  low-f reguency  controllers, 
now  be  performed  by  very  small  digital  circuits. 


can 


Figure  7.1  ADAPTI7E  COMTBOLLEB 
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Figure  7.2  DIGITAL  BLOCK  DIAGBAH 


In  converting  an  analog  controller  to  a  digital 
controller,  the  process  can  he  broken  down  into  the 
following  steps: 

1.  Determine  the  desired  analog  transfer  function. 

2.  Set  the  sampling  freguency. 

3.  Apply  the  bilinear  z-transf ormation. 

4.  Hatch  one  point  in  the  s  domain  to  the  z  domain. 

5.  Obtain  the  optimum  constant  coefficients. 

6.  Obtain  the  digital  transfer  function. 

7.  Obtain  the  simulation  diagram. 

The  optimal  controller  parameters  can  be  stored  in 
lenory.  Intel  company  markets  a  4  megabit  non-volatile 
read/ write  bubble  memory.  It  is  supported  by  a  VSLI 


controller  which  provides  a  black  box  interface.  It  is  easy 
to  use  and  can  be  used  with  any  8-  or  16-bit  aicroproces- 
sors.  Tke  bubble  memory  advantage  is: 

1.  Fast  access  tiae  compared  with  disk  or  tape. 

2.  Non-volatile. 

3.  Hide  temperature  range  of  operation. 

4.  forking  storage. 

5.  Portable  operation 

6.  Low  power. 

7-  High  reliability. 


VIII.  COHCLOSIOB5  AND  BECO EMENDATIONS  FOB  POT OB E  STOP! 

A.  CONCLUSIONS 

In  designing  the  controller,  three  different  ship  ncdels 
*ere  used.  Using  the  second-order  Nomoto  model  Equation  2.1 
allowed  comparison  of  results  with  Reid’s  [Bef.  7]  [Bef.  10] 
work.  It  is  clear  that  the  answers  obtained  by  function 
minimization  agree  closely  with  Reid’s  results  as  shown  in 
Tables  4  and  5.  A  tetter  description  of  the  ship  is  the 
third-order  Nomoto  model  which  involves  both  the  sway  and 
yaw  equations.  This  model  includes  the  two  dominating  poles 
of  the  ship.  The  best  model  to  describe  the  dynamics  of  the 
ship  is  a  Taylor’s  series  expansion.  This  allowes  both 
linear  and  nonlinear  terms  in  the  equations  of  motion  to 
affect  tie  design  of  the  controller. 

To  determine  which  controller  structure  would  provide 
the  minimum  cost  due  to  steering,  various  structures  were 
studied.  It  was  found  that  the  dynamics  of  the  plant  deter¬ 
mines  the  optimum  structure  for  the  controller.  In  calm 
water  study,  when  using  a  second-order  Nomoto  model,  the 
best  structure  was  controller  A.  When  the  third-order  Ncmcto 
model  Equation  2.2  was  used  the  best  structure  was 
controller  C,  but  the  difference  is  slight.  Observe  that  in 
each  case  the  controller  zeros  cancel  the  plant  poles,  when 
the  equations  of  motion  were  used  for  the  plant,  the  best 
structure  was  controller  C.  When  the  equations  of  motion 
were  coupled  to  a  sea  state  generator  and  the  cost  function 
was  minimized  in  the  presence  of  a  seaway,  the  best  struc¬ 
ture  was  controller  A.  This  study  concludes  that  controller 
A  should  be  used. 

A  function  minimization  subroutine  is  an  engineer's  tool 
which  can  be  used  in  many  engineering  problems.  Previously  a 


aatchcd  filter  was  designed  for  the  Naval  Postgraduate 
School  research  project  on  the  Space  Transportation  System 
(STS)  for  the  Get  Away  Special  Program.  It  was  matched  to 
the  signature  of  the  auxiliary  power  unit  (APU)  on  board  the 
space  shuttle.  The  goal  was  to  turn  on  the  solid  state 
recording  system  before  lift  off,  to  record  the  acoustic 
power  generated  inside  the  shuttle  bay.  Basically  the 
Batched  filter  is  a  finite  Impulse  Besponse  (FIE)  filter 
with  the  weights  calculated  to  obtain  the  least  squared 
error  of  the  desired  output  when  the  input  is  the  signature 
of  the  APO.  Figure  8.1  shows  the  scheme  used  to  evaluate 
the  FIS  weights. 

B.  EECCHHENDATIONS  FCE  FUTURE  STUD! 

In  the  future  most  ships  both  military  and  commercial 
will  have  GPS  receivers  as  part  of  their  navigation  equip¬ 
ment.  Using  extremely  accurate  three-dimensional  position 
and  velocity  information  from  satellite  platforms  will  allow 
ships  to  navigate  accurately  in  and  out  of  ports.  The  func¬ 
tion  linimization  subroutine  is  a  powerful  tool  for 
designing  the  controller.  This  routine  simply  takes  the 
inputs  that  require  linimization  and  adjusts  the  parameters 
to  accomplish  this  task.  The  cost  function  for  the  added 
drag  due  to  steering  is  a  function  of  yaw  error  and  rudder 
motion.  The  use  of  function  minimization  and  N AVSTAB/GPS 
provides  the  means  for  optimization  for  guidance  and  cont- 
roll.  There  are  several  areas  that  need  future  study  and 
work. 

1.  Should  the  objective  change  to  track  following  then 
it  is  necessary  to  minimize  the  yaw  error  only.  This  would 
be  very  important  both  militarily  and  commercially  should  a 
port  be  lined.  If  the  ship  could  follow  a  stringent  route, 
knowledge  of  mine  locations  would  allow  access. 


Finite  Impulse  Response  Filter 
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FUNCTION  MINIMIZATION 
SUBROUTINE 

J  =  COST  FUNCTION  =  Je2(  n  )dt 


Figure  8.1  BATCHED  FILTER  DESIGH 

2.  A  controller  for  orbit  keeping  for  satellites  with 
High-Energy  Laser  weapons  would  be  very  important.  The  snail 
far- field  spot  size  cf  a  focused  laser  bean  can  be  selec¬ 
tively  focused  on  the  aost  vulnerable  coaponent  on  the 
target.  Facilitating  precision  energy  deposition,  and 
greatly  increasing  the  probability  of  a  kill. 

3.  An  adaptive  controller  to  nininize  track  error  on 
board  a  cruise  aissile  could  be  programmed  for  selective 
targets. 


APPENDIX  A 

PROGBAH  10  CALCULATE  OPTIHAL  GAINS 

The  program  is  set  ap  to  calculate  the  optimal  gains  for 
controller  A.  It  is  referenced  in  Chapter  five  and  six.  It 
can  easily  be  modified  to  obtain  optimal  gains  for  the  rest 
of  the  controllers.  After  obtaining  the  optimal  gains  the 
program  most  be  modified  to  do  a  simulation.  The  progras  has 
sufficient  comments  for  appropriate  changes.  It  is  refer¬ 
enced  in  Chapter  tvo. 

This  program  can  be  modified  to  obtain  the  Ncmoto 
models.  It  is  referenced  in  Chapter  two.  The  following  need 
to  be  changed. 

C  GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K=XX  (1 J 
TP1=XX  (2) 

TZ=XX  (3) 

TP2=XX  (4) 

C  ESBCE  SIGNAL  TO  EEIVE  RUDDEB  (YAW  ACTUAL  -  YAW  CCMM AN I) 
c  FOB  EQUATIONS  OF  MOTION. 

E= YAW  -  YAWC 

C  ERROB  SIGNAL  TO  EEIVE  RUDDER  (YAW  COMMAND  -  YAW  ACTUAI) 

C  FCB  NCMOTO  3  ED  OEEEB  MODEL. 

D2=YAWC-YAW2 
11=  (D2-X2)  /T  P 1 
X3=K*(TZ*XUX2) 

X4=  (X3-X5)  /T  E2 
C  INTEGRATION 

I2=X2*X1*DELT 
X5=X5+X4*DEL1 
YAW2=YAW2*X5 *£ELT 
C  CCST  FUNCTION 

TDIFF=TDIFF ♦  (YAW-YAW2) **2 


PBOGHAH  TO  CALCULATE  OPTIMAL  GAINS  FOB  CONTBOILEB 


XUfI 


//GABCIA  JOB  (2220,0356) ,» RESEARCH ', CLASS=J 
//  EXEC  FRTXCLGP,IMSI=DP,BEGION=1024K 
//FCBT.SYSIN  DD  * 

C 

C  THIS  PBOGBAM  WILL  OBTAIN  TEE  CONTBOILEB  OPTIMAL  GAINS. 

C  IT  IS  BEFEBENCED  IN  CHAPTER  5. 

C 

C  IN  OBDEB  TO  PEBFOBM  SIMOLATION  ONLY  WHEN  GAINS  B AYE  BEEN 
C  OBTAINED  CHANGE  XS(*1  TO  X(*)  AND  DELETE  XU(*),AND  XL(*). 
DIMENSION  XS  (3)  ,  XU  (3)  ,  XL  (3) 

XS  (1)  =.  427 
XS  (2  =48.66 
XS  (3  =10.7 

C  XS(I)  IS  THE  STABTING  GOESS 

C  XL(I)  IS  THE  LOHEB  LIMIT  FOR  THE  I'TH  VARIABLE 
C  XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I'TH  VARIABLE 
XL  (1)  =.  1 
XU  1  =10. 

XL  2  =1. 

XU  2  =200. 

XL  3  =.10 
XU  3  =100. 

C  A  DESCRIPTION  OF  THE  FOLLOSING  PABAMETERS 
C  IS  DISCUSSED  IN  ECXPLX 
B=9./13. 

NTA=1000 

NPB=100 

NAV=0 

NV=3 

IP=0 

C  TEE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C  CAIL  PLANT  (X) 

C  IE  ONLY  SIMULATION  IS  WANTED 

CAIL  BOXPLX(NV,NAVf NPB ,NTA, B, XS, IP, XU,XL, YMN, IEB) 
WBITE  (6,25) 

25  FORMAT J1X,'  CETIMAL  GAINS',/) 

DO  30  i»1*3 


DO  30 


U\j  JV  A*  I  I J 

WRITE  (6. 40)  I  ,XS  (I) 
FOBMAT  (IX, 'X  (',12,')=', 


STOP 

END 

SUBROUTINE  PLANT (XX) 
SUBROUTINE  PLANT  (XX)  SIM 
COMMON  TDIFF 


F 1 4.  7) 


ULATES  THE  SHIP 


BEAL*8  L,L2,L3,I4,L5,L6 

REAL*8  X, XDOT,  i. YDOT, U, UDOT, V. VDOT, YA8,R,RDOT 
BEAL*8  TIME.EIIME.XUDOT.XUU.XVE.XVKxDD 
BEAL*8  YV , YB  ,  YD, Y V VR, Y VRB , YVVV , YBBB , YDDD, YVDOT 
BEAL* 8  NV.NB .NC.NVVB. NV BB, NV VV, NBBB , NDDD, NRDOT 
BE AL*8  BH0.I2,FX,FY, MZ, XP, MASS, DELT 
BEAL*8  DYAWE, z AW E, YAW C, IS E, ISB,LAMDA,D 
REAL*8  Kl,T1,T2,D,X2,Dx2,CH(11),S 
DIMENSION  XX (3) 


CLOSE  LOOP  ANALYSIS  WITH  FILTEB 

INITIAL  CONDITIONS  FOB  INTEGRATION 

SIMULATION  END  TIME  IN  SECONDS 
ETIME=600. 

TIME=0. 0 
ICCUNT= 1 

INITIALIZE  THE  COST  FUNCTION 
I SE=0. 0 
I  SR=0.  0 
TDIFF=0. 0 
LAMDA=4.2 


c 

c 

c 


c 

c 

c 

c 


ABE  FIXED  COORDINATES  ON  EABTH 


C 

C 


GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K1=XX(1) 

T1=XX(2 
T2=IX(3 
X,IDC1,Y,YD0T 
X=0.0 
Y=0.  0 
XDOT=0.  0 
YDOT=0. 0 

U,UDCT, V, VDOT  ABE  FIXED  COOBDINATES  ON  SHIP 

v=o'o 

ODOT=0. 0 
¥DOT=0. 0 
YAI=0.0 
B-0. 0 
BC01=0. 0 

CSDEEED  SPEED  IN  FEET/SEC 
54.01  FT/SEC=32  KNOTS 
UC=54.01 

AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  S PEED 
o=oc 

ANGLE 


(UC) 


ISEA=0  (CALM  BATES)  ISEA=1  (SEA  STATE) 

PABAMETEBS 


D  =  BUDDER 
D=0-  0 
1*880.5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
15=1*14 
I6=L*L5 

SEA  DISTURBANCE 

FCBCES  IN  X.Y  DIRECTION  COMPUTED  IN  POUND  FORCE 
MOMENTS  IN  2 
FX=0. 

FY=0. 

MZ=0. 

IS EA  IS  A  SWITCH 
ISEA=1 

HYDBODYNAHIC  COEFFICIENTS  ABE  INSERTED  HERE  AS 
BHO=1.9876 
MASS=  (.  0044)  *  (.5*RHO*L3) 

IZ=  (0.00028)  *  ,5*RHO*L5) 

YAWE=0. 0 
X2=0.  0 
DX2=0.0 

200  CONTINUE 

S=DSQBTj[u**2*¥**2) 

INPUT  YAH  COMMAND 
YAHC=0. 0 

IF  (TIME. GE. 0.0)  YAWC=0.0 

EBBOB  Signal  to  ebiye  rudder  (yah  actual  -  yaw  ordered) 

(  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAHC 
DX2=  (YAWE-X2)  /T2 
D=R1*]t1*DX2VX2) 

AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

XODOT  IS  THE  ADDED  MASS  TEBM  WHICH  MUST  BE  CHANGED  FOR 
DIFFERENT  ENCOUNTER  ANGLES  ,  SPEED  ,  ENCOUNTER  FREQUENCY 

XUDOT= (-.000 1 )  *  ( . 5*RHO*L3) 

XU=  (-0.  0253)  *  (.5*RHO*L2*S) 

XUU=  (-0.0003)  *  (.  5*RHO*L2) 

X¥B=  0. 0039) *  (.5*RHO*L3) 

X  VV=  -.  0012)*  ].5*RHO*L2) 

XDD=  -0.0005)  *  (.5*RHO*LZ*S**2) 

LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

YV=  (-0.  00758)  ♦(.  5*BHO*L2*S) 

IB=  0.0023)  *  (.5*RHO*L3*S) 

YD=  0.00145)  <  (,5*RHO*L2*S**2) 
l¥¥H=  (0.01)  *(.5*BHO*L3/S) 
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on  o  o  on  non  nnnnnn  n  ononn  n  nnnnn 


Y  VBB= 

Y  VVV= 
YBRR  = 
YEDD= 

YODOT  IS 


-0.003)  *  (.  5*RH0*L4/S) 

-  -  -.5*RHO*L2/S) 

5*RHO*L5/S) 

(.5*RHO*I2*S**2) 

THE  ADDtE  HASS  TEBH  WHICH  HOST  BE  CHANGED  FOB 


—  U •  UuO  I  *  .  D' 

-0.03)  *  {.5* 
0.003i  *  .5*: 

-  0.  0005  )  ♦  (.  ! 


DIFFERENT  ENCOUNTER  ANGLES,  SPEED  ,  ENCOUNTER  FREQUENCY 
Y VDOT=  {-0.0039)*(.  5*RHO*L3) 

SPEED=32  KNOTS. ENCOUNTER  ANGLE  =150  , ENCOUNTER  FREQ  =.75 
YVD0T=-2. 3043E+C6 

HCHENT  ABOUT  Z-AXIS  HYDBODYNAHIC  COEFFICIENTS 
NV=  (-0.  00213)  *  {.  5*RHO*L3*S) 

NB=  (-0.  00  105)  *  |.  5*RHO*L4*S) 


(YAW) 


NEBB=  (-0.006)  *  (•  5*RHO*L6/S) 

NDDD=?O.OOOli *  . 5*RHO*L3*S**2) 

NBEOT  IS  THE  ADDED  INERTIA  TERH  WHICH  HUST  BE  CHANGED  FOR 
DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 


NRDOT= {- 0.000 
SEEEE=32  K 


30 


27)  *(.5*RH0*L5) 

«NOTS. ENCOUNTER  ANGLE  =150  , ENCOUNTER  FREQ  =.75 
NEDOT=-1 , 5096  E  +  1 1 
SETS  SEA  STATE  TO  ZERO 

IF  (ISEA. EQ. 1)  GO  TO  30 
FX=0. 

F  Y=0. 

HZ=0. 

GC  TO  35 

TABLE  LOOK  UP  OF  SEA  DISTURBANCE. 

UNIT  12  HAS  THE  SEA  STATE  EATA  NAMED  CH 
IT  HUST  BE  SYNCHECNIZED  BY  APPROPRIATELY 
CALLING  CH  IN  THE  PROPER  TIHE  IN  THE  LOOP. 

TEE  SEA  DATA  IAS  CREATED  FOR  600  SECONDS 
WITH  AN  INCREHENTAL  INTERVAL  OF  1  SECOND. 


READ  (12)  CH 
FX=CH  '■* ' 

F Y=CH | 


35 


Si 


HZ=CH 
CCNTINUE 
U  ACTUAL  SPEED 
DC  CCHHANDED  SPEED 
XE  =  PROPELLER  THFUST 
XP=-XUU*UC**2 
EQUATIONS  OF  HOTICN 

FOR  CONSTANT  SPEED  COMMENT  THE  NEXT  TWO  INSTRUCTIONS 
UDOT=  (  (XVR  +  MASS) * V*R  ♦  XUU*U**2  ♦  XVV*V**2 
1  ♦  XDD*D*D  ♦  FX  ♦  XP  [/(MASS-XUDOT) 

VDOT= {Y V* V  ♦  (YR-HASS *□) *R  ♦  YD*D  ♦  YVVR*V«'*2*H 

1  ♦  YVHR*V*R**2  ♦  YVVV *V**3 

2  ♦  YRRR *R **3  ♦  YDDD*D**3  ♦  FY  ) / (HA SS-YVDOT) 

BEOT=l  NV*V  ♦  SR*R  ♦  ND*D  ♦  NVVR*V**2*R 

1  ♦  NVSR*V*R**2  ♦  NVVV*V**3 

2  ♦  NRRR*B**3  ♦  NDDD*D**3  ♦  HZ  )/(IZ-NRDOT) 

WEEN  TO  PRINTOUT 

IF  llCOUNT.EQ.11)  GO  TO  50 
GC  TO  300 

CCNVEBT  RADIANS  TC  DEGREES 
50  YA1DEG=  YAW*57.296 
RDEG=R*57. 296 
RDDEG=RDOT*57 .296 
DDEG=D*57. 296 
YAKC=YAWC*57«  2  96 

WRITE  (6,100)  TIHE, IP, X ,XDOT , Y, YDOT 
1  ,OC,U.UDOT, V.VDOT.YAWC, YAiDEG, RDEG, RDDEG, DDEG 
100  FCRHAx'f  lX,  '  TlHx=  *  .  Fo.  3 .  '  SEC  XP=».F10.2,'  LBF  X=  ' 
1,F8.2,*  FT  XECT=* ,F8. 4 , *  FT/SEC  Y=',F8.2, 


nnnnnnoonnnnnnonnnnnnnnnnnonnnnnnr> 


2'  PI  YDOT=*  , F8. 4, *  PT/SEC' ,/,2x. •  UC=',F8.4, 

3'  FT/SEC  U= ' ,F8. 4, *  FT/SEC  UDOT= • . Fl 0 .6 . 

4  *  FT/SEC**2  V='.F8.4,'  PT/SEC  VDOT= ' , P 10. 6, 

5»  FT/SEC* *2 ' , /,2X. ' YAHC=' , F8. 4 ,  '  DEG  YAW='  F15.7, 
6'  DEG  YAH  H ATE=* ,F1 5. 7, *  DEG/SEC  YAH  ACCE1=' 

7, FI 5. 7, *  DEG/SEC**2',/,2X, 'RUDDER  =  ',F15.7,'  DEG  • 
ICCUNT=  1 


/) 


C  TEST  IF  WANT  TO  STOP 
300  IF  (TIME.  GE.  ETIME)  GO  TO  400 
C  INTEGRATION  STEP  SIZE  DE1T 
DELT=1. 0 
C  ISIEGBATION 

0=U+UDOT*DELT 

V=V+VDOT*DELT 

R=R*RDOT*DELT 

YAW=YAfi*R*DELT 

X2=X2*DX2*DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EABTH 
XDOT=U*DCOS  (YAH)  -V*DSIN  (YAH) 

- - {yahJ 


400 


YDOT=U*DSIN  (YAH) +V*DCOS 
X=X*XDOT*DELT 
Y=Y+YDOT*DELT 
IIME-TIHE*DEL7 
ICCONT=ICOUNT*1 
ISE=ISE  ♦  LAMEA*YAWE**2 
ISB=ISH  ♦  D**2 
GO  TO  200 

J=  TDIFF  =  COST  FUNCTION 
TDIFF=ISE*ISR 


HBITE  (6.500)  I£E.ISfi,TDIFF,X 1,T1,T2 
500  POBH AT ( '  •  IX, 'I5£=  •  .FI  5.  7  ,  •  ISR=',F15.7  •  TOT  AL=1 
1gF15.7,2x,*Kl=',F15.7,2X, *T1  =  ' ,F15.  2X,  •  5:2=' ,F15. 7) 


BETUBN 
END 

DELETE  ALL  THE  FOLLOWING  SOBBOUTINE  IF  SIMOLATICN  ONLY 
AND  NOT  OPTIMIZATION  IS  WANTED 


SCEBOUTINE  BOXELX  (CATEGOBY  HO) 

PUBPOSE 

BOIPLX  IS  A  SUBROUTINE  USED  TO  SOLVE  THE  PROBLEM  CF 
locacting  A  MINIMUM  (OB  MAXIMUM)  OF  AN  AB3ITBABY  OBJECT- 
ive  function  SUBJECT  TO  ARBITRARY  EXPLICIT  AND/CR 
implicit  constraints  by  tHE  COMPLEX  METHOD  OF  M.J.  BOX. 
explicit  constraints  are  dEFINED  AS  UPPEB  AND  LOHEB 
bounds  on  the  independent  variables  IMPLICIT  constraints 
nay  be  arbitrary  function  of  the  varlABLES.  THO  FUN- 
ction  subprogram  to  evaluate  the  objective  FUNCTION  AND 
implicit  conSTRAINTS,  RESPECTIVELY.  must  be  SUPPLIED 
by  the  user  (see  EXAMPLE  BELOH) .  BOXPLX  ALSO  HAS  tHE 
option  to  perform  integer  programming,  uhere  the  values 
of  the  independent  variables  are  restricted  to  integers. 

OSAGE 

CALL  BOXPLX  (NV, N AV, NPB ,NTA, R, XS, IP, XU , XL , YMN  ,IEB) 


DESCRIPTION  OF  PARAMETEBS 

NV  AN  INTEGER  INEUT  DEFINING  THE  NUMBER  OF  INDEPENDENT 
VARIABLES  OF  THE  OBJECTIVE  FUNCTION  TO  BE  MINIMIZED. 
NOTE:  MAXIMUM  NV  +  NAV  IS  PRESENTLY  50.  NAXIMIM  NV  IS 

25.  IF  THESE  LIMITS  MUST  BE  EXCEEDED,  PUNCH  A  SOURCE 
DECK  IN  THE  USUAL  MANNER,  AND  CHANGE  THE  DIMENSION 
STATEMENTS. 


-y- 


C  HA V  AH  INTEGER  INPUT  DEFINING  THE  NUMBER  OF  AUXILIARY  var 
C  iaELES  THE  USEB  WISHES  TO  DEFINE  FOR  HIS  OWN  CONVENIENCE. 
C  TYPICALLY  HE  MAY  WISH  TO  DEFINE  THE  VALUE  OF  EACH  INELICI 
C  CONSTBAINT  FUNCTION  AS  AN  AUXILIARY  VABIABLE.  IF  THIS 
C  IS  DONE.  THE  OPTIONAL  OUTPUT  FEATURE  OF  BOXPLX  CAN  BE 
C  USED  TO  OBSERVE  THE  VALUES  OF  THOSE  CONSTRAINTS  AS  THE 
C  SCIUTIC N  PROGRESSES.  AUXILIARY  VARIABLES,  IF  USED, 

C  SHOULD  BE  EVALUATED  IN  FUNCTION  KE  (DEFINED  BELOW). 

C  NAV  HAY  BE  ZERO. 

C 

C  NPB  INPUT  INTEGEfi  CCNTBOLLING  THE  FREQUENCY  OF  ODTP.UT 
c  desired  for  diagnostic  purposes. 

C  IF  NPB  -LE.  0.  NO  CUTPOf  WILL  BE 

C  PRODUCED  BY  BOXPLX.  OTHERWISE,  THE  CURRENT  COMPLEX  OF 
C  K=  2* NV  VERTICES  AND  THEIR  CENTROID  WILL  BE  OUTPUT  AFTER 
C  EACH  NPB  PERMISSIBLE  TRIALS.  THE  NUMBER  OF  TOTAL  TRIALS, 
C  NUMBER  OF  FEASIBLE  TRIALS,  NUMBER  OF  FUNCTION  EVALUATIONS 
C  AND  NUMBER  OF  IMPLICIT  CONSTRAINT  EVALUATIONS  ARE  IN- 
C  CIUDEB  IN  THE  OUTPUT. 

C  ADDITIONALLY,  (WHEN  NPR  .GT.  0)  THE  SAME  INFORMATION 
C  WILL  BE  OUTPUT: 


IF  THE  INITIAL  POINT  IS  NOT  PEASIBLE, 

AFTER  THE  FIRST  COMPLETE  COMPLEX  IS  GENERATED, 

IF  A  FEASIBLE  VERTEX  CANNOT  BE  FOUND  AT  SOME  TRIAL, 
IF  THE  OBJECTIVE  VALUE  OF  A  VERTEX  CANNOT  BE  MADE 
NO-LCNGEE-HORST. 


t)  IF  THE  LIMIT  ON  TRIALS  (NTA)  IS  REACHED  AND, 
i)  WHEN  THE  OBJECTIVE  FUNCTION  HAS  BEEN  UNCHANGED  FOB 
>*NV  TRIALS.  INDICATING  A  LOCAL  MINIMUM  HAS  BEEN 


C  2*NV  TRIALS,  INDICATING  A  LOCAL  MINIMUM  HAS  BEEN 
C  FOUND. 


C  IF  TEE  USER  WISHES  TO  TRACE  THE  PROGRESS  OF  A  SOLUTICN, 
C  A  CHOICE  OF  NPR  =  25,  50  OR  100  IS  RECOMMENDED. 

C  NTA  INTEGER  INPUT  OF  LIMIT  ON  THE  NUMBER  OF  TRIALS 
c  allowed  in  the  calculation. 

C  IF  THE  USER  INPUTS  NTA  .LE.  0,  A  default 
C  VALUE  OF  2000  IS  USED.  WHEN  THIS  LIMIT  IS  REACHED 
C  CONTROL  RETURNS  TO  THE  CALLING  PROGRAM  WITH  THE  BEST 
C  ATTAINED  OBJECTIVE  FUNCTION  VALUE  IN  YMN,  AND  THE  BEST 
C  ATTAINED  SOLUTION  POINT  IN  XS. 


C  R  A  REAL  NUMBER  INPUT  TO  DEFINE  THE  FIRST  RANDOM  NOMEER 
C  USED  IN  DEVELOPING  THE  INITIAL  COMPLEX  OF  2*NV  VEfiTICIES. 
C  lO.  .GT.  R  .LT.  1.)  IF  R  IS  NOT  WITHIN  THESE  BOUNDS, 

C  IT  HILL  BE  REPLACE!  BY  1./3.  . 


C  IT 
C 

C  XS 


C  XS  INPUT  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV+NAV. 
c  the  first  nv  oust  contain  a 
C  FEASIBLE  ORIGIN  FCB  STARTING  THE  CAL- 

C  CULATICN.  THE  LAST  NAV  NEED  NOT  BE  INITIALIZED.  UPON 
C  HETUBN  FROM  BOXPLX,  THE  FIBST  NV  ELEMENTS  OF  THE  ARRAY 
C  CONTAIN  THE  COORDINATES  OF  THE  MINIMUM  OBJECTIVE 
C  function.  AND  THE  REMAINING  NAV  (NAV  .GE.  0)  CONTAIN  THe 
C  values  of  THE  CORRESPONDING  AUXILIARY  VARIABLES. 

C 

C  JP  INTEGER  INPUT  FOR  OPTIONAL  INTEGER  PROGRAMMING. 

C  if  ip=1,  THE  VALOIS  OF  THE  INDEPENDENT  VARIABLES  WILL 
C  fc€  replaced  WITH  INTEGER  VALUES  (STILL  STORED  AS  REALM)  . 


C 

C  XU 


c 

C  XL 


A  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUTTING  THE 
er  BOUND  ON  EACH  INDEPENDENT  VARIABLE.  (EACH  EXPLICIT 
STBAINT).  INPUT  VALUES  ARE  SLIGHTLY  ALTERED  BY  BOXPLX. 


3RED  BY  BOXPLX. 


C  XL  A  BEAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUTTING  THE 
c  lower  bound  on  each  independent 
C  VABIABLE,  (EACH  EXPLICIT  CONstraint)  . 


V  *  -v  N  *.  NS.*.  \ 


I 
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BOTE:  FOR  BOTH  XU  AND  XL  CHOOSE  REASONABLE 

VALUES  IF  NONE  ABE  GIVEN,  NOT  VALUES  WHICH  ARE 
aagaitudes  ABOVE  CB  BELOV  THE  EXTECTED  SOLUTION, 
infut  values  are  SLIGHTLY  ALTERED  BY  BOXPLX. 


YHB  THIS  OUTPUT  IS  THE  VALUE 
fuacTICN,CORRESPO  BUNG  TO  THE 


(REAL*4) 

SOLUTION 


OF  THE  OBJECTIVE 
POINT  OUTPUT  IN  XS 


IEB  INTEGER  ERROR  RETURN.  TO  BE  INTERROGATED  UPCN 
return  FROM  BOXPLX.  IER  HILL  BE  ONE  OF  THE  FOLLOWING: 

=-1  CANNOT  FIND  FEASIBLE  VERTEX  OR  FEASIBLE  CENTROID 
AT  THE  START  OR  A  RESTART  (SEE  *  METHOD’  BELOV). 

=0  FUNCTION  VALUE  UNCHANGED  FOR  • N'  TRIALS.  (WHERE 
N=6*NV* 1 0)  THIS  IS  THE  NORMAL  RETURN  PARAMETER. 

=1  CANNOT  DEVELCE  FEASIBLE  VERTEX. 

=2  CANNOT  DEVELCE  A  NO-LCNGEB-WORST  VERTEX. 

=3  LIMIT  ON  TRIALS  REACHED.  (NTA  EXCEEDED) 

NOTE:  VALID  RESUl'IS  MAY  BE  RETURNED  IN  ANY  OF  THE 

ABOVE  CASES. 

EXAMPLE  OF  USAGE 


THIS  EXAMPLE  NINIMIZ 
the  EXTERNAL  FUNCTIC 
varlAELES  X(1)  S  X  (2 
function  X(3)  S  X(4) 
variables  (see  EXTER 


E  MINIMIZES  THE  OBJECTIVE  FUNCTION  SHOWN  IN 
L  FUNCTION  FE(X).  THERE  ARE  TWO  INDEPENDENT 
(1)  &  1(2),  AND  TWO  IMPLICIT  CONSTRAINT 
3)  S  X  (4)  WHICH  ARE  EVALUATED  AS  AUXILIARY 
see  EXTERNAL  FUNCTION  KE(X)  ). 


DIMENSION  XS  (4)  #X0  (2)  ,XL  (2) 


STARTING  GUESS 
XS(1)  =  1.0 
XS J2)  =  0.5 
UPPER  LIMITS 
XU  (1)  =  6.0 
XU  (2)  =  6.0 
LOWER  LIMITS 
XI(1)  =  0.0 
XL  (2)  =  0.0 


R  =  9./13. 
NTA  =  5000 
NPR  =  50 
NAV  =  2 
NV  =  2 
IP  *  0 


CALL  BOXPLX  (NV, NAV, NPR .NTA, R, XS.IP, XU,XL,YMN,IEE) 
WRITE76, 1 )  (  (XS (I)  .1=  1,4)  , YHN, IER) 

1  FORMAT  (///),'  THE  POINT  IS  LOCATED  AT  (XS(I)=)  • 


1  FORMAT  (////,'  THE  POINT  IS  LOCATED  A' 
2. 4 (e 13. 7,5x) ,//, 

3*  AND  THE  FUNCTION  VALUE  IS  ',E13.7,» 


cc 

c 

STOP 

c 

END 

cc 

cc 

c 

FUNCTION 

FUNCTION  KE (X) 

EVALUATE  CONSTRAINTS.  SET  KE=0 
is  viCLATED,  OR  SET  KE=1  IF  ANY 
constraint  is  violated. 
DIMENSION  X (4) 


IF  NO  IMPLICIT 
IMPLICIT 


-  *.15) 


CONSTRAINT 


DIMENSION 
XI  =  X ( 1 ) 
X2  =  X  (2) 
KE  =  0 
XJ3)  =  XI 
?*...<*  <3 )  . 


X  (4) 


♦  1.732051*X2 
,LT.  _0.  .OR.  X  (3)  .GT. 


1/1.732051  -X2 


GO  TO 
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IF  (1(4)  .  GE.  C.)  RETURN 

1  KI  =  1 
RETURN 
END 


C  DIHENSION  X (4) 

cc 

CC  THIS  IS  THE  OBJECTIVE  FUNCTION. 

C  *(9--(X(1)-3.)**2)/(46. 76538)) 

C  BETURN 

C  END 

C 

C  METHOD 

C 

C  THE  COMPLEX  METHOD  IS  AN  EXTENSION  AND  ADAPTION  OF 
c  the  simple  method  cf  linear  programming. 

C  STARTING  WITH  ANY  CNE  feasible  point  in  n-dimension 
C  A  "COMPLEX"  OF  2*N  vertices  is  constructed  by 
C  SELECTING  RANDOM  ECINTS  WITHIN  THE  feasible 
C  BEGIC3.  FOR  THIS  PURPOSE  N  COORDINATES  ARE  FIRST 
C  RANDOMLY  CHOSEN  WITHIN  THE  SPACE  BOUNDED  BY  EXPLICIT  CCN- 
C  STRAINTS.  THIS  DEFINES  A  TRIAL  INITIAL  VERTEX, 
c  it  is  then  checked  for  possible  violation 
C  OF  IMPLICIT  CONSTRAINTS.  IF  one  or  more  are  violated, 

C  THE  TRIAL  INITIAL  VERTEX  IS  DISPLACED  half  of  its 
C  DISTANCE  FROM  THE  CENTROID  OF  PREVIOUSLY  SELECTED  initial 
C  VERTICES.  IF  NECESSARY  THIS  DISPLACEMENT  PROCESS  IS 
C  REPEATED  UNTIL  THE  VERTEX  HAS  BECOME  FEASIBLE.  IF  THIS 
c  FAIL  TO  happen  after  5*n*10  displacements, 

C  THE  SOLUTION  IS  AEANDONED.  AFTER  EACH  VERTEX  IS  ADDED 
C  TO  THE  COMPLEX,  THE  CURRENT  centroid  is  checked  fcr 
C  FEASIBILITY.  IF  IT  IS  INFEASIBLE,  the  last  trail 
C  VERTEX  IS  ABANDONED  AND  AN  EFFORT  TO  GENERATE  an  alter- 
C  ATIVE  TRIAL  VEHTEX  IS  MADE.  IF  5*N*10  VERTICES  ARE 
C  ABANDCNED  CONSECUTIVELY,  THE  SOLUTION  IS  TERMINATED. 

C 

C  IF  AN  INITIAL  COMPIEX  IS  ESTABLISHED,  THE  BASIC 
c  computation  loop  is  initiated. 

C  THESE  INSTRUCTIONS  FIND  THE  CURRENT  WORST  vertex,  that 
C  IS,  THE  VERTEX  WITH  THE  LARGEST  CORRESPONDING  value  for 
C  THE  OBJECTIVE  FUNCTION,  AND  REPLACE  THAT  VERTEX  BY 
C  ITS  OVEB-REFLECTICN  THROUGH  THE  CENTROID  OF  ALL  OTHER 
c  vertices,  (if  the  vertex  to  be 
C  REPLACED  IS  CONSIDERED  AS  A  VECTOR  IN  n-space, 

C  ITS  CVEB-REFLECTICfi  IS  OPPOSITE  IN  DIRECTION,  IN- 
C  CR EASED  IN  LENGTH  EY  THE  FACTOR  1.3,  AND  COLLINEAB  WITH 
C  THE  REPLACED  VERTEX  AND  CENTROID  OF  ALL  OTHER  VERTICES.) 

C 

C  WHIN  AN  OVER— REFLECTION  IS  NOT  FEASIBLE  OR  REMAINS 
c  WORST.  IT  IS  CONSIBEBED  NOT-PERMISSI BLE 
C  AND  IS  DISPLACED  EALFWAY  TOWARD  THE  CENTROID. 

C  AFTER  FOUR  SUCH  ATTEMPTS  ARE  MADE  UNSUCCESSFUL! X 
C  EVERY  FIFTH  ATTEMPT  IS  MADE  BY  REFLECTING  THE  OFFENDING 
C  VERTEX  THROUGH  THE  PRESENT  BEST 

C  VERTEX,  INSTEAD  OF  THROUGH  THE  CENTROID.  IF  5*n+10 
C  DISPLACEMENTS  AND  OVER-REFLECTIONS  OCCUR  WITHOUT  A 
C  SUCCESSFUL  (PERMISSIBLE)  RESULT,  THE  CURRENT  BEST 
C  VERTEX  IS  TAKEN  AS  AN  INITIAL  FEASIBLE  POINT  FOR  A 
C  RESTABT  RUN  OF  THE  COMPLETE  PROCESS. 

C  RESTARTING  IS  ALSO  UNDERTAKEN  WHEN  6*nv*10  CONSECUTIVE 
C  TRIALS  HAVE  BEEN  MADE  WITH  NO  SIGNIFICANT  CHANGE  IN  THE 
C  VALUE  OF  THE  OBJECTIVE  FUNCTION.  IN  ALL  CASES, 

C  RESTARTING  IS  INHIBITED  IF  THE  LAST  RESTART  DID  NOT 
C  PRODUCE  A  SIGNIFICANT  IMPROVEMENT  IN  THE  MINIMUM 
c  ATTAINED. 
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C  IT  IS  RECOMMENDED  THAT  THE  OSEB  BEAD  THE  BEFEBENCE  POB 
C  FUBTHEB  USEFUL  INFORMATION.  IT  SHOULD  BE  NOTED  THAT  THE 
C  ALGOBITHH  DEFINED  THEBE  HAS  BEEN  ALTEEED  TO  FIND  THE 
C  CONSTRAINED  MINIMUM,  BATHER  THAN  THE  MAXIMUM. 

C 

C 

C 

C  REMARKS 

C 

C  THE  INTEGER  PROGRAMMING  OPTION  HAS  ADDED  TO  THIS  EROGRAM 
C  AS  SUGGESTED  IN  BEFEBENCE  (2).  A  MIXED 
c  integer/continuous  variable  version  of  boxplx 
C  MOULD  BE  EASY  TO  CBEATE  BY  DEclaring  "ip"  to  be  an  array 
C  OF  NV  CONTROL  VARIABLES  WHERE  IP  (if = 1  would  indicate 
C  THAT  THE  I-TH  VARIABLE  IS  TO  BE  CONFINED  to  integer 
C  VALUES-  EACH  STATEMENT  OF  THE  FORM  •IF  (IP  .  EQ. T)  '  etc. 

C  WOULD  THEN  NEED  TO  BE  ALTEEED  TO  'IF  (IP  (I)  .  EQ.  1)'  etc. 
C  ,  WHEBE  THE  SUBSCEIPT  IS  APPROPRIATELY  CHOSEN.  NORMALLY, 
C  XU  AND  XL  VALUES  ABE  ALTERED  TO  BE  AN  EPSILON  'WITHIN' 
c  actual  values 

C  DECLARED  BY  THE  USER.  THIS  ADJUSTMENT  IS  NOT  HADE 
C  WHEN  IP=1. 

C 

C  NOTE:  NO  NON-LINEAR  PBOGR AMMING  ALGORITHM  CAN  GUARANTEE 
c  that  the  answer  found  is  the  global 
C  MINIMUM,  RATHER  THAN  JUST  A  local  minimum,  however, 

C  ACCORDING  TO  REF. 2,  THE  COMPLEX  method  has  an  advantage 
C  IN  THAT  IT  TENDS  TO  FIND  THE  GLOBAL  minimum  more 
C  FREQUENTLY  THAN  MANY  OTHER  NCN-LI NEAR  PROGRAM- 
C  MING  ALGORITHMS. 

C 

C  IT  SHOULD  BE  NOTED  THAT  THE  AUXILIARY  VARIABLE  FEATURE 
C  CAN  ALSO  BE  USED  TO  DEAL  8ITH 

C  PRCELEMS  CONTAINING  EQUALITY  CONstraints.  any  equality 
C  CONSTRAINT  IMPLIES  THAT  A  GIVEN  VARiable  is  not  truly 
C  INDEPENDENT.  THEREFORE,  IN  GENERAL,  ONE  variable 
C  INVOLVED  IN  AN  EQUALITY  CONSTRAINT  CAN  BE  RENUMBERED  from 
C  THE  SET  OF  NV  INDEPENDENT  VARIABLES  AND  ADDED  TO  THE  SET 
C  OF-  HAV  AUXILIARY  VARIABLES.  THIS  USUALLY  INVOLVES 
C  renumbering  THE  INDEPENDENT  VARIABLES  OF  THE  GIVEN 
C  problem 

c  Subroutines  and  functions  required 

C  SUBROUTINE  'BOUT'  AND  FUNCTION  ' FBV'  ARE  INTEGRAL 
C  parts  of  THE  BOXPIX  PACKAGE. 

C 

C  TWO  FUNCTIONS  MUST  BE  SUPPLIED  BY  THE  USER.  THE  FIRST, 
c  ke  (xl,  is  used  to  evaluate  the  implicit 

C  CONSTRAINTS.  SET  KE=0  AT  THE  beginning  of  the  function 
C  .  TEEN  EVALUATE  THE  IMPLICIT  CONstraints.  in  the  example 
C  AECVE,  THE  FIRST  CONSTRAINT,  X  (3 },  must  be  within  the 
C  RANGE  (0.  .LE.  X  ( 3)  -LE.  6.).  THE  SECOND  constraint  x  (4) 
C  ,  MUST  BE  .GE.  0.  .  IF  EITHER  CONSTRAINT  IS  not  within 
C  THESE  BOUNDS.  CONTEOL  IS  TRANSFERRED  TO  STATEMENT  1, 

C  AND  KE  IS  SET  TO  "1"  AND  CONTROL  IS  RETURNED  TO  BOXPIX. 

C 

C  THE  SECOND  FUNCTION  THE  USER  MUST  PROVIDE  EVALUATES  TIE 
c  objective  function,  it  is 

C  CALLED  FE  (X)  AS  SHOWN  IN  THE  EXAM  pie  above,  and  fe 
C  MOST  EE  SET  TO  THE  VALUE  OF  THE  OBJECTIVE  function 
C  CORRESPONDING  TO  CURRENT  VALUES  OF  THE  NV  INDEPENDENT 
C  VABIAELES  IN  ARRAY  'X*. 

C 

C  REFERENCES 

C 
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SUBROUTINE  BOXPLX  |NV,NAV,NPR/NTZ/RZ/XS,IP,BU,BL/YMN/IER) 

DIMENSION  V  (50,50),  FUN  (50)  ,  SUM(2S),  CEN  (25)  ,  XS(NV) 
1 , lu  {nv)  ,  bl  (NV) 

KV  =  5 
EP  =  1.  E-6 
NTA  =  2000 

IF  (NTZ.GT.0)  bTA  =  NTZ 
H  =  RZ 

IF  (R.LE.O..OR.R.GE.  1.)  R=1./3. 

NVT  =  NV+NAV 

TOTAL  VARS,  EXPLICIT  PLUS  IMPLICIT 

NT  *  0 

CURRENT  TRIAL  NO. 

NET  =  0 

CURRENT  NO.  OF  PERMISSIBLE  TRIALS 
NTFS  =  0 

CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGEE 
CHECK  FEASIBILITY  OF  START  POINT 


DC  4  1=1, NV 
VT  =  XS(I) 

IF  _(BL|f)  .LE.  VT)  GO  TO  1 

VT  =  BL  (I) 

GC  TO  2 

IF  (BU  (I)  .  GE  .  V  T)  GO  TO  3 
II  =  I 
VT  =  BU  (I) 

IF  (NPR.GT.0)  WRITE  {6.45 


2  IF  (NPR.GT. 

3  V  IX.  II  =  VT 


0)  WRITE  (6,49)  II 


VJ  1,1)  =  VT 

CEN  (I)  =  VT 

IF  jlP. EQ.  1)  GC  TO  4 

Bl  (I)  =  BL(I)  +AMAX1  (EP, 

BU  ( ii  =  BUll) -AMAX1  (EP, 

SUM  (I)  =  VT 


EP* ABS 
EP  * ABS 


NCE  =  1 

NUMBER  CF  CONSTRAINT  EVALUATIONS 

1=1 

IF  (KE(V  (1,  1)  1.EQ.0)  GO  TO  5 


IF  (KE(V  (1,  1)  ).EQ.0)  GC 
IF  JNPR.LE.O)  GO  TO  12 
WRITE  (6,50) 


WRITE  (6 
GC  TO  12 
5  NFE  =  1 


RUKEER  OF  VERTICES  (K)  =  2  TIMES  NC.  OF  VARIABLES. 
K  =  2*N V 

NUMBER  OF  DISPLACEMENTS  ALLOWED. 


NLIfl  =  5*NV+ 1 0 


NUMBEB  OF  CONSECUTIVE  TRIALS  KITH  UNCHANGED  FE  TC 
ter linate. 

NCT  =  NLIM+N7 
ALPHA  =  1.3 
FF  =  K 
FKM  =  FK-1. 

BETA  =  ALPHA+  1. 

INSURE  SEED  OF  RANDOM  NUHEEB  GENERATOR  IS  ODD. 

XCH  =  R  *  1  •  E7 

IF  (MOD  (IQR,2).EQ.O)  ICE=IQB+101 

SET  UP  INITIAL  VEBTICES 
FUN  ( 1)  =  FE(V(1,1)) 

„  YBN  =  FUN  (1) 

6  FI  =  1. 

FUNCLD  =  FUN  ( 1) 

DC  15  1=2, K 
FI  =  FI  + 1  . 

LIMT  =  0 

7  LIMT  =  LIMT+1 

END  CALCULATION  IF  FEASIBLE  CENTROID  CANNOT  BE  FOUND 
IE  (LIMT. GE. NLIM)  GO  TO  11 


DC  8  J=1,NV 


£ ANDCM  NUMBER  GEKEBATOS  (BANDU) 

ICR  =  IQR*65539 

IF  (IQR.LT.O)  IQB  =  IQH  +  2  147483647+ 1 
BOX  =  IQR 

BOX  =  RQX*.  4 6  566 1 3 E- 9 

VjJ,I)  =  BL  (J)+RQX*(BU(J)-BL(J)) 

c  ▼(J+D-AliT'fFtJ/lJV.S) 

8  CONTINUE 


i.  =  flU*¥.H03u0  I  Jt“3 

r,I)  =  BL  (J)+RQX*(BU(U)-BL(J)) 
(it.  EQ.  1)  V  (J,I)  =AINT  (V  (j/lj  V.  5) 
ITINUE 


DO  10  L= 1 ,  NLI M 
NCE  =  NCE+1 

IF  (KE(V(1,I)  )  .EQ.O)  GO  TO  13 
DO  9  J=  1 ,  N  V 

VI  =  .5*  (V(J,I)+CEN(J)) 

IF  (IP.EQ.1)  VI  =  AINT  (VT  +  .5) 


IF  (IP.EQ.1)  VI  =  AINT 
V  (J.I)  =  VT 
S  CONTINUE 

10  CONTINUE 

11  IF  (NPR.LE.O)  GO  TO  12 


IF  (NPR.LE.O)  GO  TO  12 
BBITE  (6,51)  I 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 


12  IER  =  -1 
GO  TO  48 

13  DO  14  J= 1 ,N V 


•  W  It  KJ  —  I  m  a  1 

SUM  (J)  =  SUM  (J)+V (J,I) 
14  CEN(J)  =  SUM(0)/FI 


TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 
NCE  =  NCE+1 


IF  (KE(CEN)  .EC.O)  GO  TO  60 
SUMjJ)7=  SUM  ( J)  -7(J,I) 


60  NEE  =  NFE+1 

FUN  (I)  =  FE  ( V  (1,1)) 
15  CONTINUE  " 


K 
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END  CF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

IF  (NPR.LE.O)  GO  TO  17 

CALX  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,G) 

FIND  THE  WORST  VERTEX,  TEE  *J«TH. 

J  =  1 

DC  16  I - 2  .  X 

IF  (FUN  (J)  .GE.EDN  (I)  )  GOTO  16 

16  CONTINUE 

EASIC  LOOP.  ELIMINATE  EACH  WORST  VERTEX  IN  TURN, 
it  Bust  become  NC  LONGER  WORST.  NOT  MERELY  IMPROVED, 
find  next-to- verte x,  THE  'JN'TH  ONE. 

17  JN  =  1 

IF  (J.EQ.  1)  JN  =  2 
DC  18  1=1. K 


Lf  \J  I  VJ  J.  —  I  /  XV 

IF  (I.EQ.J)  GC  TO  18 

IF  (FUN  (JN)  -GE.FUN  (I)  )  GO  TO  18 


JN  =  I 
18  CONTINUE 

II MT  =  NUMBER  OF  MOVES  DURING  THIS  TRIAL  TOWARD  THE 
centroid  DUE  TO  FUNCTION  VALUE. 

LIMT  =  1 

CCMEUTE  CENTROII  AND  OVER  REFLECT  WORST  VERTEX. 


DC  19  1=1  .NV 
VI  =  V(I,J) 

SUM  (I)  =  SUM  { I)  -VT 

CEN(l[  =  SUM(l[/FKM 

VT  =  BETA*CEN  (I) -ALPHA*VT 

IF  (IP.EQ.1)  \I  =  AINT  (VT+.5) 

INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AMAX1  (A  MINI  (VT,BU  (I)  )  ,  BL  (I)  ) 

NT  =  NT+1 

CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 

20  DC  25  N=1,NLIH 
NCE  =  NCE ♦ 1 

IF  (KE  (V  (1,  J)  )  .EQ.O)  GO  TO  26 

EVERY  ’ KV ' TH  TIME,  OVER-EEFLECT  THE  OFFENDING  VERTEX 
through  the  BEST  VERTEX. 

IF  t«OD  (N,KV)  .  KE.  0)  GO  TO  22 
CA1L  FBV  (K,FU6,M) 

DC  21  1=1. NV 

VT  =  BETA*V  (I  . B)  - ALPH A*V  (I,  J) 

IF  (IP.EQ.1)  vl  =  AINT  (VT  +  .5 

21  V(I,J)  =  AMAX  1  (AMIN  1  ( VT , BU  (I)  )  ,BL  (I)) 

GC  TO  24 

CONSTRAINT  VIOLATION:  MCVE  NEW  POINT  TOWARD  CENTROID. 

22  DC  23  1=1, NV 


yw  X.  J  4.  -  1  f  H  V 

VI  =  .5*  JCEN  (I)*V  (I,  J)  ) 

IF  (IP.EQ.1)  VI  =  AINT  ( VT  ♦.  ! 


Vjl.J)  =  VT 
23  CONTINUE 
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nn  n  noon  nn  nn 


24  NT  =  NT+1 

25  CONTINUE 
C 

IIR  =  1 

C  CANNOT  GET  FEASIELE  VERTEX  BY  MOVING  TOR ARD  CENTROID, 

C  CB  EY  OVER-REEL  ICTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE.O)  GO  TO  42 
WRITE  (6,52)  NT, J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GC  TO  42 

FEASIELE  VERTEX  FOUND,  EVALUATE  THE  OBJECTIVE  FUNCTION 

26  NFE  =  NFE  +  1 
FUNTRY  =  FE  (V  (1,  J)  ) 

TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED. 

AFO  =  ABS (FUNTEY-FUNOLD) 

ASX  =  AMAxI  (AES  (£P*FU  NOID)  ,  EP) 


ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS  FOR  DIAGNOSTIC 
purposes  only. 

WRITE  (6,99)  J, AFO, ASX, FUNTRY, FUNOLD, FUN (J) , FUN (JN) 
1,NTFS,N 

99  FORMAT  ( 1 X, 13 ,6E 1 5. 7 , 21 5) 

IF  (AFO.GT.AMX)  GO  TO  27 
NTFS  =  NTFS+1 
IF  (NTFS.LT.  NCT)  GO  TO  28 
IER  =  0 

IF  (NPR.LF.O)  GO  TO  42 
WRITE  { 6,53}  K 

CALI  BOUT  (NT  ,NPT, NFE,NCE,NV, NVT,V,K,FUN,CEN,0) 

GO  TO  42 
27  NTFS  =  0 


IS  TEE  NEW  VERTEX  NO  LONGER  WORST? 

28  IF  (FUNTRY.LT. FUN  (JN)  )  GO  TO  34 

THIAL  VERTEX  IS  STILL  WORST:  ADJUST  TOWARD  CENTROID. 
EVERY  'KV'TH  TICE,  OVER-REFLeCT  THE  OFFENDING  VERTEX 
through  the  BEST  VERTEX. 

LIST  =  LIMT+1 


IF  (HOD  (LIST,  KV)  .NE.O)  GO  TO  30 
CALI  FBV  (K,FD*,M) 


DC  29  1=1, NV 

VT  =  BETA*V  (I -M) -ALPHA*  VJI,J) 

IF  (IP.EQ.1)  VT  =  AINT  ( VT  +  .  5  j 

29  V(I,J)  =  AM  AX  1  (A  MINI  (VT,BU  (I)  )  , 

GC  TO  32 

30  DO  31  1=1, NV 
VT  =  .5*  (CEN  (I)+V  (I,  J)  ) 

IF  (IP.EQ.1)  \1  =  AINT  ( VT  +  .5) 

VJI.J)  =  VT 

31  CONTINUE 


BL  (I)) 


32  IF  (LIST. IT. NLIH)  GO  TO  33 

C  CANNOT  MAKE  THE  »J'TH  VERTEX  NO  LONGER  WORST  EY 
c  displacing  toward 

C  THE  CENTROID  OR  BY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 
IER  =  2 

IF  (NPR  .LE.  0)  GO  TO  42 
WRITE  (6,52)  NT.  J 

CALI  BOUT  (NT , NPT ,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 

GC  TO  42 

33  NT  =  NT+' 


on  nnn  on  onnonnn  on  non  nnoonnoo  non  on 


GO  TO  20 

SUCCESS:  V  E  HAVE  A  REPLACEMENT  FOE  7  EHT2X  J. 

34  FUN  (J)  =  FUNTRX 
FUNOLD  =  FUNTEY 
NPT  =  NPT+1 

I7ERY  100 1 TH  PERMISSIBLE  TRIAL.  RECOMPUTE  CENTROID 
suanation  to  AVCID  CREEPING  ERROR. 

IF  (MOD  (NPT,  ICC)  .NE.  0)  GO  TO  37 

DO  36  1=1. NV 
SUM  (I)  =  0. 

DC  35  N  =  1  ,  K 

35  SUM  (I)  =  SUM  { I)  +  V  (I,  N) 

CEN(I)  =  SUM  (I) /FK 

36  CONTINUE 

LC  =  0 
GO  TO  39 

37  DO  38  1=1, NV 

38  SUM  (I)  =  SUM  (I)+V  (I,  J) 

LC  =  J 

39  IF  (NPR.LE.0)  GO  TO  40 

IF  (MOD  (NPT,NFB) -NE.  0)  GO  TO  40 

CALL  BOUT  (NT , NPT, NFE , NCE, NV,NVT,V,K,FUN,CEN,LC) 

EAS  THE  MAX.  NUMEER  OF  TRIALS  BEEN  REACHED  WITHOUT 
convergence?  iF  NOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 

NEXT-TO-WOBST  VERTEX  NOW  BECOMES  WORST. 

J  =  JN 
GO  TO  17 

41  IEB  =  3 

IF  (NPR.GT.0)  SHITE  (6,54) 

COLIECTOR  POINT  FOR  ALL  ENDINGS. 

1)  CANNOT  DEVELCE  FEASIBLE  VERTEX.  IER  =  1 

2  CANNOT  DEVELOP  A  NO- LONGER-WORST  VERTEX.  IER  =  2 

3  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  IER  =  0 

4;  LIMIT  ON  TRIALS  REACHED.  IER  =  3 

5  CANNOT  FIND  FEASIBLE  VERTEX  AT  START.  IER  =  -1 

42  CONTINUE 

FIND  BEST  VERTEX. 

CALL  FBV  (K  .FUN,  M) 

IF  (IER.GE.3)  GO  TO  44 

RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  THAN 
the  previous,  OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  (NPR.LE.0)  GO  TO  43 
WRITE  (6,55)  (F.YMN,  FUN  (M)  1 

43  IF  (FUN  (8)  .GE.YMN)  GO  TO  47 

IF  (ABS  (FUN (M) -YMN) . LE. AMAX1 (EP, EP*YMN)  )  GO  TO  47 

GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

44  YMN  =  FUN  (M) 

FUN  ( 1)  =  FUN  (M) 

DO  45  1=1, NV 
CENfl)  =  V  (I ,  N) 

SUM  (I)  =  7(1, B) 
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45  V  (1,1)  =  V(I,K) 

DC  46  1=1, NVT 

46  XS(I)  =  V  (I,  M) 


IP 
47  IF 


IF  (IER.LT.3)  GO  TO  6 
IF  jNPR.LE.O).  GO  TO  48 

CALI  BOOT  (NT,NPT,NFE,NCE,NV,NVT,V, K,FUN,V  (1,M) ,-1) 
WRITE  (6,56)  FUN  (A) 


48  RETURN 

49  FORMAT  (50H0INDEX  AND  DIBECTION  OF  OUTLYING 
1  variable  at  start i5) 

50  FORMAT  (50  HO  IMPLICIT  CONSTRAINT  VIOLATED  AT 
Istart.  dead  end.) 

51  FCBMAT  ( ' OC  AN NCT  FIND  FE ASIBLE ' , 14, • TB  VERTEX  OB 
Icentroid  at  start.*) 

52  FCBMAT  (10H0AI  TRIAL  I4,54H  CANNOT  FIND  FEASIBLE 
Ivertex  which  is  no 

2ICNGEB  WORST.I4, 15X,  '  RESTART  FROM  BEST  VERTEX.') 

53  FORMAT  (4 OHOFUNCTION  HAS  BEEN  ALMOST  UNCHANGED 
Ifor  i5,7h  trails) 

54  FCBMAT  (27H0LIMIT  ON  TRIALS  EXCEEDED.  ) 

55  FORMAT ( * OBEST  VERTEX  IS  NO.1, 13, 'OLD  MIN  WAS',E15.7, 
1  '  NEW  MIN  IS  ' ,E15. 7) 

56  FORMAT  ('OMIN  OBJECTIVE  FUNCTION  IS  »,E15.7) 

END 

SUBROUTINE  FBV  (K,FUN,M) 

DIMENSION  FUN  (50) 

H  =  1 

DO  1  1=2, K 

IF_  (FUN  (M)  .LE.  EUN  (I)  )  GOTO  1 

1  CCNTINUE 

BETUBN 

END 

SUBROUTINE  BOUT  (NT,NPT ,BFE, NCE.NV, NVT,V,K,FN,C, IK) 
DIMENSION  V(50,5<3),  FN(50),  C(25) 

WHITE  (6,4)  N1,NPT, NFE, NCE 

DC  1  1=1,  K 

WRITE  (6.5)  FN(I)  ,  (V(J,I)  ,J=1,NV) 

IF  (HVT.LE.HV)  GO  TO  1 
N VP  =  NV+1 

WRITE  (6,6)  (V  (J,I)  ,J=NV£,NVT) 

1  CONTINUE 

IF  (IK.NE.O)  GO  TO  2 

<C(I),I=1,NV) 

BET  URN 

2  IF  (IK.GE.O)  GC  TO  3 
WRITE  (6,8)  (C  (I)  ,  1=  1  ,  N V) 

RETURN 

RETURN 

4  FORMAT  ('ONO.  TOTAL  TRIALS  =  ',I5,4X, 

1'no.  feasible  trails  =  ',i5,4x, 

2 ' NO.  FUNCTION  EVALUATIONS  =  *,I5,4X, 

3'nc.  constraint  evaluations  =  *,r5/ 

4*0  FUNCTION  V ALU E ' , 6X , ' INDEPEND ENT  VARIABLES/ 

5dependENT  OB  IMPLICIT  CONSTRAINTS') 

5  FORMAT  ( 1 H  , E 1 8.7 . 2X, 7E 14 . 7/  (2 IX, 7E 1 4. 7) ) 

6  FORMAT  2 IX, 7E  14.7) 

7  FCBMAT  10H0CENTROID  1 1 X, 7  El  4 .7/ (2 1 X. 7E1 4. 7) ) 

8  FCBMAT  '0  BEST  VERT  EX • , 7X, 7E1 4 . 7/ (2  IX, 7 E 14. 7) ) 


-  -S  . 


9^FCRHAT  (* OCENTEOID  LESS  7 X • , 12 ,2X, 7E 14.7/ ( 2 1 X, 
END  * 

FUNCTION  FE(X) 

DIMENSION  Xj3) 

CCHMCN  TDIFF 
CALL  PLANT (X) 

FE=TDIFF 

BETUHN 

END 

FUNCTION  KE (X) 

DIMENSION  X  (3) 

KE=0 

EETUEN 

END 

//GO.SYSIN  DD  * 

/* 

//GO. FT 12F00 1  DD  DISf=SHR, DSN-MSS. S2160. A34 1 


APPENDII  B 

EXAHPIE  PBOBLEB  OSIHG  ICSOS 


The  purpose  of  this  example  is  to  demonstrate  the 
performance  of  the  program.  Consider  the  control  system  of 
Figure  4.1  with  controller  C.  Figure  B. 1  shows  the  block 
diagram  to  evaluate  the  controller  parameters. 


CONTROLLER  I  NOMOTO  THIRO  ORDER  MODEL 


Ihe  differential  equations  describing  the  system  and  its 
desired  performance  aie: 

2*2=  DX2 
I4=BX4 
X5=DX5 
B  =DB 
YAJJ=B 

Defining  the  following  cost  function: 

J=  y*(LAMDA*yAWf**2+D**2)  dt 
Defining  the  special  functions: 

YAHE=YAWC-YAW 
D12=  (YAHE-X2)  /12 
X3=K1 * (T 1 *DX2  +X2) 

D  X4=  (X3-X4)  /T4 
d=  (t3*dx4+x4) 
dx5=  (d-x5)/tp1 
x8=k*  (tz*dx5  +  x5) 
di=  (x8-r)  /tp2 
Defining  the  constants: 

YAHC=1.0 
K=. 14771 
1Z= 11.2833 
1F1=6.4699 
TE2=53.  7931 
LANCA=4.2 

and  using  YAiC=1.0  the  optimal  solution  found  by 
the  program  is: 

K1=. 4179916 
11=53.69932 
12=4.970023 
13=6.294369 
14=13.85919 
C0S1=68. 04735 
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Table  34  shows  the  specifications  of  this  problea  with 
the  free  parameter  cptiaua  values  found.  Figure  B.2  shows 
the  actual  yaw  and  rudder  response. 


TABLE  34 
ICSOS  OUTPUT 


- SPECIFICATIONS- 


£  INITIAL  CCNCIT IONS: 


VARIABLES 
X2  *  .0 
X4  *  .0 
X5  *  ,0 
R  =  .0 
YAW  =  .C 
J  =  .0 
TIME  =  .0 


CONSTANTS: 

YAw£  *  l.OOOOOCOOO 
K  =  .147710000C 
TZ  =  11.283300CO 
TP1  =  6.469900000 
TP2  =  53*  793 100C0 
LAMGA  *  4.20COCCOOQ 

FREE  PARAMETERS: 

K1  I  CV=  .4179916  LL= 

T1  :  oy*  53.69928  LL* 

T2  :  0 V=  4.970C29  L L  = 
T3  :  0V=  6.294369  LL* 

T4  :  0V=  13.  85917  LL* 

SPECIAL  FUNCTICNS: 

YAWE  *  YAWC-YAV* 

DX2  *  (YAWE-X2J/T2 
X3  *  Kl*  (T1*DX2+X2  ) 

DX4  *  IX3-X4J/T4 
D  =  (T3*CX4+X4J 
0X5  *  ( C-X5) /TF1 
X8  *  K*(TZ*D  X5+X5) 

DR  *  (X8-RI/TP2 

DERIVATIVES: 

D(X2  /D (TIME  )  *  * 

0X2 

0CX4  /□ ( TIME  )  *  = 

0X4  x 

D(X5  /D (T IME  I  =  = 

0X5 

D(R  /D ( T IME  )  =  = 

OR 

DIYAW  /C( TIME  »  *  * 

R 

DCJ  /0( TIME  )  *  * 

LAM0A*YAWE**2+D**2 


.1CCC000 

1.000000 

l.OOCOOO 

l.OOCOOO 

l.CCCOOO 


LL* 
LL* 
UL® 
UL  = 
LL* 


1.00000  0 
ICO.  000  0 
!C.  0000  0 
10.00000 
20. 0000 C 


OUTPUTS: 

TITLE:  ACTUAL  YAW  ANC  RUDDER  RESFCNSE 
TABULATE:  T  IFF.  D  R  YAW 

AT  INTERVAL  2.CCOCOOOOO 
PLOT:  0  YAW 

AGAINST:  TIME  AT  INTERVAL  2.000000C00 
END  CALCULATION  WHEN  TIME  .GE.  600.000 


AfA  ■♦♦♦♦* 


Figure  B.2  YAH  AMD  BODDEB  VS.  TIME 
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